DUD ;HY KNOX LIBRARY 

NAV; I,. POSTGRADUATE SCHOOL 

MONTEREY CA 93943-5101 




Approved for public release; distribution is unlimited. 



Anode Sheath Contributions in Plasma Thrusters 

by 

John F^Riggs 

Lieutenant Commander, United States Navy 
B.A., University of Kansas, 1982 

Submitted in partial fulfillment 
of the requirements for the degrees of 

AERONAUTICAL & ASTRONAUTICAL ENGINEER 

and 

MASTER OF SCIENCE IN ASTRONAUTICAL ENGINEERING 

from the 

NAVAL POSTGRADUATE SCHOOL 
March 1994 



1 

REPORT DOCUMENTATION PAGE 


Form Approved 0MB No. 0704 

_____ 


Public repotting burden for this collection of information is estimated to average 1 hour per response, including the time for revie\Mng instruction, 
searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments 
regarding this burden esiin.ata or other aspect of coikcikm of infof ii.alioTi, iadudi**g suggestions for reducing ibis burden, to Washiugon 

headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202^302, and 
to the Office of Management and Budget, Paperwork Reduction Project (0704-0188) Washington DC 20503. 


1. AGENCY USE ONLY (Leave fc/an/tj 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 

March 1994 Engineer’s and Master’s Thesis 


4. TITLE AND SUBTITLE ANODE SHEATH CONTRIBUTIONS IN 
PLASMA THRUSTERS 


5. FUNDING NUMBERS 


6. AUTHOR(S) Riggs, John Forrest 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 
Naval Postgraduate School 
Monterey CA 93943-5000 


8. PERFORMING 
ORGANIZATION 
REPORT NUMBER 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDR£SS(ES) 


10. 

SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


11. SUPPLEMENTARY NOTES The views expressed in this thesis are those of the author and do not 
reflect the official policy or position of the Department of Defense or the U.S. Government. 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 
Approved for public release; distribution is unlimited. 


12b. DISTRIBUTION CODE 
A 


13. ABSTRACT {maximum 200 words) 



Contributions of the anode to Magnetoplasmadynamic (MPD) thruster performance are considered. 
High energy losses at this electrode, surface erosion, and sheath/ionization effects must be controlled 
in designs of practical interest. Current constriction or spotting at the anode, evolving into localized 
surface damage and considerable throat erosion, is shown to be related to the electron temperature's 
(TJ rise above the gas temperature (T„). An elementary one-dimensional description of a collisional 
sheath which highlights the role of is presented. Computations to model the one-dimensional sheath 
' are attempted using a set of five coupled first-order, nonlinear differential equations describing the 
electric field, as well as the species current and number densities. For a large temperature 
nonequilibrium (i.e., > > TJ, the one-dimensional approach fails to give reasonable answers and 

a multidimensional description is deemed necessary. Thus, anode spotting may be precipitated by the 
elevation of among other factors. A review of transpiration cooling as a means of recouping some 
! anode power is included. Active anode cooling via transpiration cooling would result in (1) quenching 
Tj, (2) adding "hot" propellant to exhaust, and (3) reducing the local electron Hall parameter. 
However, significant technical problems remain. 

15. 

NUMBER OF 
PAGES 106 

16. 



PRICE CODE 





18. 


19. 


20. 


SECURITY CLASSIFI- 


SECURITY CLASSIFI- 


SECURITY CLASSIFI- 


LIMITATION OF 


CATION OF REPORT 


CATION OF THIS PAGE 


CATION OF ABSTRACT 


ABSTRACT 


Unclassified 


Unclassified 


Unclassified 


UL 



NSN 7540-01-280-5500 



1 



Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std. 239-18 



ABSTRACT 



Contributions of the anode to Magnetoplasmadynamic (MPD) thruster performance 
are considered. High energy losses at this electrode, surface erosion, and sheath/ionization 
effects must be controlled in designs of practical interest. Current constriction or spotting 
at the anode, evolving into localized surface damage and considerable throat erosion, is 
shown to be related to the electron temperature’s (TJ rise above the gas temperature (TJ. 
An elementary one-dimensional description of a coliisional sheath which highlights the role 
of Tj is presented. Computations to model the one-dimensional sheath are attempted using 
a set of five coupled first-order, nonlinear differential equations describing the electric field, 
as well as the species current and number densities. For a large temperature nonequdibrium 
(i.e.,T,> >T„), the one-dimensional approach fails to give reasonable answers and a 
multidimensional description is deemed necesssary. Thus, anode spotting may be 
precipitated by the elevation of among other factors. A review of transpiration cooling 
as a means of recouping some anode power is included. Active anode cooling via 
transpiration cooling would result in (1) quenching T., (2) adding "hot" propellant to exhaust, 
and (3) reducing the local electron Hall parameter. However, significant technical problems 



remain. 






TABLE OF CONTENTS 

I. INTRODUCTION 

II. LITERATURE REVIEW 

III. ANODE DESCRIPTION 

A. THRUSTER GEOMETRY DESCRIPTION 

B. ELEMENTARY SHEATH FORMULAE DESCRIPTION 

1. Discussion 

2. Simplified Formulation 

3. Approximate Formulation 

a. Effects of Temperature on Anode Constriction 

(1) Case I: T, = T, = T„ (Equilibrium) 

(2) Case II: T, > > T, = T„ (Two-Temperature) . . . . 

b. Similarity to Vacuum Arc Phenomena 

C. COMPUTER CODE 

D. COMPUTATIONAL RESULTS 

E. ANODE FALL VOLTAGE 

IV. TRANSPIRATION COOLING 

V. CONCLUSIONS AND RECOMMENDATIONS 



1 

5 

7 

7 

11 

11 

14 

21 

22 

23 

25 

27 

29 

30 

36 

38 

43 



IV 



DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-5101 



APPENDIX - APPLICABLE FORTRAN PROGRAMS 46 

LIST OF REFERENCES 92 

INITIAL DISTRIBUTION LIST 96 



V 



LIST OF TABLES 



I. NOMENCLATURE 15 

LIST OF FIGURES 

1. Magnetoplasmadynamic (MPD) Thruster 8 

2. Space Plasma Thruster 10 

3. Electric Field Between Two Electrodes 11 

4. Electric Field and First Two Space Derivatives 19 

5. Electric Field and Species Approximation 24 

6. Two-dimensional Model of Current Paths 26 

7. Anode Discharge Modes 28 

8. Ionization Coefficient v as a Function of E/N 30 

9. Species Number Density Plots for Individual Computer Run 32 

10. Electric Field as a Function of y 34 

11. Species Number Density as a Function of y 35 



VI 



ACKNOWLEDGMENT 



I dedicate this work to tny dear wife Lin, for her support and understanding 
during my absence as a geographical bachelor of almost three years. 

My sincere thanks to my parents for encouraging me to pursue my ambitions from 
childhood on, and to Oscar Biblarz, Ph.D. , for patiently helping me to understand plasma 
physics. 



Vll 



I. INTRODUCTION 



Several types of space flight propulsion systems have been developed over the 
years. These include chemical, nuclear, electric and solar propulsion. The majority 
of space thrusters to date have been chemical rockets. Although the Chinese used 
rockets over 800 years ago, true development of rocket propulsion took place during 
this century [Ref. 1]. Chemical thrusters give high thrust-to-weight ratios, larger than 
unity, and have been fully developed in the form of space launch vehicles and 
attitude control thrusters. In contrast, other propulsion systems have been developed 
only to the proof-of-concept stage, and essentially remain at this stage of 
development. Nuclear propulsion was studied with the NERVA (Nuclear Engine for 
Rocket Vehicle Application) thruster in the 1960’s, and abandoned [Ref. 2:pp. 517- 
519]. Electric propulsion flights during the 1960’s included the U.S. SERT-1 (Space 
Electric Rocket Test) in 1964 and the U.S.S.R. Yantari-1 rocket in 1966. Solar- 
electric propulsion was demonstrated via the SERT-2 rocket in 1970, powering the 
electric thruster from power generated by solar cells. Further electric propulsion 
research flights in the 1980’s included the U.S. Navy’s NOVA-1 satellite in 1981, and 
Japan’s MS-T4 satellite,launched from the Space Shuttle. Beyond this, nonchemical 
thrusters have only been used in auxiliary roles, such as station-keeping and attitude 
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control on geosynchronous satellites. NASA’s Project PATHFINDER in the mid- 
1980’s proposed the use of a megawatt-level electric plasma thruster for a manned 
Mars mission. However, development of this project was never funded. 

In comparing the different propulsion schemes, a primary performance indicator 
is specific impulse, defined as the ratio of thrust to the rate of propellant usage, or 
alternately, propellant effective exhaust velocity (u^), divided by the sea-level 
gravitational constant, (go). 



I = ^ 



u 

— sec 

§o 



( 1 ) 



Chemical rockets are inherently limited in performance by the total energy available 
in the fuel/oxidizer combustion process, so that the total enthalpy available for 
conversion into exhaust kinetic energy is limited. Exhaust velocity is also limited by 
material heating limitations of the combustion chamber and nozzle throat, and 
"frozen flow Losses" (unrecoverable energy deposition in internal modes of the gas) 
[Ref. 3:pp. 4-5]. Peak specific impulse for liquid chemical propellants is presently on 
the order of 450 seconds. This capability is completely sufficient for the tasks of 
launch to low earth orbit (LEO), attitude control, station keeping, and such missions. 
However, for missions such as manned interplanetary exploration, chemical 
propulsion can be shown to be clearly inadequate. A comparative analysis of a Mars 
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mission using chemical and electric propulsion systems shows the large difference in 
mass payload ratio (final mass/initial launch mass) for the two systems. A chemical 
system using a high impulse Hohman ellipse trajectory delivers a maximum of 
approximately 10% to 18% of the launch mass to the Martian surface [Ref. 4:p. 115]. 
In comparison, an electric system using a low impulse spiral trajectory could deliver 
from 20% to 60% of the launch mass, depending on the desired transit time. Each 
mission assumes transit from low Earth orbit to Mars orbit. An electric propulsion 
system would still need a high thrust propulsion system to reach the Martian surface 
Ref. 5:pp. 344-346]. The large difference in payload ratio is due to the much larger 
exhaust velocity and more efficient use of fuel by electric propulsion. Thus, some 
form of electric or hybrid electric thruster would seem to be in order for such 
interplanetary missions. However, due to the low thrust-to-weight ratio of electric 
thrusters, they must be launched into orbit by other means. Their usefulness is 
restricted to space thrusters, not to launch systems. 

With specific impulses of as high as 10,000 seconds, electric propulsion offers 
the performance envelope needed for manned interplanetary missions. Electric 
propulsion is divided into three types of thrusters: electrostatic, electrothermal, and 
electromagnetic. The type relevant to this work is the magnetoplasmadynamic (MPD) 
thruster, an electromagnetic propulsion system that utilizes the Lorentz force created 
by an electric current together with its induced magnetic field to propel a gas that 
has been heated to the plasma state. According to electromagnetic theory, a 
conductor carrying a current produces an induced magnetic force perpendicular to 
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the current. The applied electric field and its induced magnetic field interact to 
produce the Lorentz force (F = JxB) perpendicular to both fields on the 
conductors. This briefly summarizes the concept behind the "self-field" MPD 
accelerator [Ref. 2:pp. 485-486]. MPD performance is enhanced by adding magnetic 
coils to the thruster, thus strengthening the magnetic field and, as a consequence, the 
Lorentz force and thrust. This thruster is appropriately called an "applied-field" MPD 
thruster. MPD thrusters have shown specific impulses of up to 7,000 seconds and 
efficiencies as high as 70% [Ref. 6:pp. 2-3]. Performance of MPD thrusters is limited 
by several factors, including electrode erosion, current spotting, frozen flow losses, 
and electrode power deposition. Specifically, anode power deposition is the single 
largest power loss mechanism in MPD thrusters operating at submegawatt power 
levels [Ref. 7]. In the following work, we review and analytically model the MPD 
anode, including the sheath and anode potential drop. 
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II. LITERATURE REVIEW 



Anode losses significantly limit magnetoplasmadynamic (MPD) thruster 
performance. Much effort has been placed on characterizing these losses and on the 
nature of power deposition in the anode [Refs. 8-14]. As much as 80% of thruster 
total power may end up being deposited in the anode at sub-megawatt power levels 
[Refs. 8,15]. This power deposition together with current constriction at the anode 
surface present serious problems to thruster cooling and performance, as well as to 
anode lifetime. Before any practical design can be achieved, a more thorough 
understanding of the phenomena at the anode, particularly the anode sheath, must 
be gained. Studies have shown that the anode power fraction depends on thruster 
power, current, mass flow rate, and the parameter J^/ih [Refs. 8,12,13,16]. It has also 
been shown that the anode fall voltage is inversely proportional to anode current 
density [Refs. 13,16]. It is believed that a better understanding of the role of an 
elevated electron temperature, of current flow dimensionality, and of current 
unsteadiness are prerequisites for the evolution of any practical MPD thruster design. 

Computer codes that accurately describe observed data from steady-state MPD 
thrusters have been developed [Refs. 17-19]. However, these codes do not adequately 
describe observed data from quasi-steady thruster experiments. It has been suggested 
that the lack of proper electrode modelling (i.e., sheaths and fall potentials) in these 
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codes may explain this discrepancy [Ref. 6]. Limited analytical work has been done 
in modelling the sheath and ambipolar regions at the anode, influenced perhaps by 
the difficult set of coupled, nonlinear partial differential equations involved. Hugel 
[Ref. 12] and Subramaniam [Ref. 20] address the influence of the sheath region, but 
do not model the electric field, temperature, or sheath fall voltage. 

Given the minuscule extent of the sheath versus thruster anode curvature, the 
problem at first appears one-dimensional in nature. A one-dimensional, collisional, 
equilibrium solution can satisfactorily reproduce the observed electric field and 
charge density distributions for the entire sheath and ambipolar regions for a sheath 
where the electron temperature equals that of the heavy species [Ref. 21]. However, 
this model cannot describe any decrease in current density away from the surface, or 
current constriction, at the anode surface which might be necessary in 
nonequilibrium. A two-dimensional model, developed by Biblarz and Dolson [Ref. 
14], represents these phenomena and predicts the voltage drop in the region. It is 
shown that the sheath must account for a majority of the anode voltage drop, and 
that the sheath extent must be greater than the Debye length [Refs. 14,21]. Thus, 
a combination of one- and two-dimensional approaches appears to better describe 
sheath behavior. Incorporation of modelling of this sort may improve the ability of 
the computer codes cited above to properly describe quasi-steady thrusters. 

Next, a description of the anode region is presented in order to delineate some 
of the possible effects of temperature. 
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III. ANODE DESCRIPTION 



A. THRUSTER GEOMETRY DESCRIPTION 

The majority of plasma thrusters to date have consisted of a central cathode 
rod surrounded by an annular shell anode, as shown in Figure 1 [Ref. 23]. The 
thruster illustrated is sufficient to produce needed thrust at current levels above one 
kiloamp. Below this level, an external magnetic field produced by an annular magnet 
is needed to ensure sufficient Lorentz force on the plasma propellant to meet thrust 
requirements. [Ref. 8]. As illustrated in Figure 1, the JxB body force simplifies into 
an axial ( body force, which provides direct electromagnetic thrust ("blowing"), 
and a radial ( -JzBq) body force, which provides electromagnetic compression of the 
plasma and a subsequent pressure force along the cathode surface ("pumping"). [Ref. 
6 ] 

A notable exception to this geometry is the Stationary Plasma Thruster (SPT), 
a design from the former Soviet Union. The SPT is an example of a plasma 
propulsion system known as a Hall Current Plasma accelerator. An electric field is 
applied axially to a stream of flowing plasma, in addition to a magnetic field with a 
strong radial component, which is applied by an external electromagnet. When the 
axial electric field is applied and a current flows through the plasma, an azimuthal 
component of current is induced, i.e., the "Hall" current. 
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Figure 1 - Magnetoplasmadynamic (MPD) Thruster, with Axial and Radial Forces 
on Plasma Indicated. (Ref. 23] 
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Thrust is produced by electrostatically accelerating the ions in the plasma, as well as 
through the induced Lorentz force mentioned above. A strong radial magnetic field 
is applied to the plasma, whose properties are controlled to make the electron 
Larmor radius’ small compared to the mean free path^ while the ion Larmor radius 
is comparatively large. As a consequence, electron mobility in the axial direction is 
greatly reduced. Thus, the electric field energy is given mainly to the ions, producing 
axial ion acceleration. Collisions with neutral particles serve to accelerate the entire 
neutral plasma. [Ref. 24] 

A pair of the final prototype design developed, the SPT-100, have been 
acquired by NASA recently from Fakel Enterprise in Kaliningrad, Russia, and are 
undergoing performance evaluation at the Jet Propulsion Laboratory. Designed at 
the Kurchatov Institute of Atomic Energy (lAE) in Moscow, USSR in the 1960’s, 
smaller versions of the SPT-100 (SPT-50 and SPT-70) were flown beginning in 1972\ 
A specific impulse of 1,600 seconds and 50% efficiency, as well as space flights of 
fifty similar thrusters is claimed. The specific operational characteristics of the 
thruster are not well understood presently. Bohm diffusion of electrons and a 
phenomenon called "near-wall conductivity" have been proposed to explain the 
thruster’s operation. This thruster is shown in Figure 2. [Ref. 25] 

‘ Larmor radius is the radius of the helix traversed by a charged particle moving 
in a magnetic field. 

^ Mean free path is the distance traveled by a particle before making a collision. 

^The suffbc (i.e.,"-70") indicates the characteristic diameter of the thruster in 
millimeters. 
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Figure 2 - Stationary Plasma Thruster [Ref. 25] 
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B. ELEMENTARY SHEATH FORMULAE DESCRIPTION 



1. Discussion 

Voltage losses and anode power deposition account for most of the 
inefficiency of plasma thrusters. In order to understand these losses, the anode region 
must be understood and related phenomena explained and modelled. As shown in 
Figure 3, a substantial drop in voltage occurs in a short distance from the anode 
surface. 







-Positive Column--^ 



Cathode 



Anode 



Figure 3 - Electric Field Between Two Electrodes, Including "Positive Column". [Ref. 
27] 

The anode fall region may be divided into two parts, the sheath and ambipolar 
regions. The plasma attempts to adjust itself near electrodes so as to shield the main 
body of the plasma from the electric field [Ref. 26]. The sheath is the region closest 
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to the anode surface within which the ion and electron number densities are unequal, 
with the electrons dominating the region. A high electron space charge exists at the 
anode surface. This is caused by the anode collecting incoming electrons in 
completing the arc current with the cathode. Positive ions are produced within the 
sheath by electron impact of neutral gas molecules, and the ions are repelled toward 
the cathode. At the cathode end of the anode drop region, the density of positive 
ions is high enough to almost neutralize the electron space charge, thus forming the 
positive column or core plasma. The essential positive ion current is created in this 
way near the anode. A more complete description of this process may be found in 
Cobine (Ref. 27] and von Engel [Ref. 28]. A fundamental characteristic of plasma 
behavior is its tendency toward electrical neutrality. Whenever local charge 
concentrations arise or external potentials are introduced into a system, these are 
shielded out in a distance known as the "Debye length"\ This distance must be much 
smaller than the system dimension for the ionized gas to be considered a plasma 
[Ref. 29]. Equation (2) gives the Debye length [Ref. 26]. 






N 



e kT 

= 69.0 

ne^ 






T 

- (m) 

n 



( 2 ) 



The Debye length effectively describes the radius of a shell around a charged 
particle outside of which the potential of the particle is not seen. 
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Another distance of interest is the electron mean free path, or distance traveled by 
a particle before making a collision. Equation (3) is from a derivation of Lin, Resler, 
and Kantrowitz [Refs. 30,31] giving the mean free path, with being the 
approximate sheath length. 



X 



0.12 



/ 



V 



1 

n^(eV3kT)Mn(A.eV3kT) 



(3) 



Since the sheath extends at most a few mean-free lengths from the anode surface, 
curvature of the anode does not affect the governing equations in high pressure 
discharges. Thus, the region may be described in one dimension, the distance "y" 
from the anode surface. While the Debye length is sometimes assumed as the sheath 
extent. Reference 22 showed that the sheath thickness is a function of the anode fall 
voltage and the electron temperature. Equation (4) gives the appropriate form. 



X 



S 





(4) 



An example case with a fall voltage of 100 volts gives a sheath extent of 
X^ = 2.352x10'^ m. This compares to a computed Debye length of 
X^ = 1.69 0x10"® m. Therefore, the sheath can be an order of magnitude larger 
than the Debye length. 
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Nasser [Ref. 32] discusses an elementary theoretical approach to the glow 
discharge problem. He suggests a set of four one-dimensional ordinary differential 
equations, including the electron and ion current and number density equations, in 
addition to Poisson’s equation. Most solution attempts have failed, with the boundary 
conditions being identified as the culprit. A similar attempt for the plasma thruster 
is discussed below. 

2. Simplified Formulation 

The steady probe equations are first written [Ref. 21] in their simplest 
form. The anode is assumed to operate as a heavily biased probe, which is true for 
low enough currents when the anode is not a source of ions. Whenever the 
temperature can be considered fbced, the energy equations are implicitly satisfied 
and, since ion inertia is neglected, the resulting set consists only of two species 
continuity equations and Poisson’s equation. These equations are written in terms 
of y, which is the coordinate outward from the planar positive surface. Constants and 
variables are listed in Table 1. 
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TABLE 1 - NOMENCLATURE 



a... characteristic length of plasma 


n„. ..species number density at core plasma 


Di e... species diffusion coefficient 


N... total number density 


e...elementary charge constant 


T... temperature 


E...electric field 


Tq... neutral species temperature 


Eq ...electric field at anode surface 


a ...two-body recombination coefficient 


E„...electric field at core plasma 


Eg. ..permittivity constant 


ji g. ..species current density 


V ...ionization coefficient 


J... total current 


Pi e— species mobility coefficient 


k...Boltzmann’s constant 


anode fall potential 


K...current parameter 


X... mean-free distance 


Hi g...species number density 


X^... Debye length 


iie...time rate-of-change of rig 


Xg... Sheath thickness 



Note: Species subscripts denote ions (i) and electrons (e). 
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dn 

I = e/Li nE - (eD) — ^ 
J. , V ./ 


(5) 


dn 

j. = e^inE + (eD^)^ 

dy 


(6) 


dE e / . 

— = — (n - n ) 

dy e/ ' 


(7) 


J = j, 


(8) 


eD 

11 = 

kT 

i.e 


(9) 



Here, the j’s are species contributions to the total current density . The existence of 
negative charges as free electrons is pivotal in the formulation. Next, the Einstein 
relation, equation (9), is introduced to write the mobilities in terms of the diffusion 
coefficients. We assume that the diffusion coefficients remain constant in the 
problem. 

Equations (5) and (6) are next solved for dn^ Jdy. The species current 
density equations are found from the net reaction rate of the plasma. Equations (10) 
and (11) combine to produce space derivatives for species current density. 
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n = V n - an n 

c I c 1 < 



(10) 




= eri 



(11) 



dy dy ' 

Combining equations (5)-(ll) produces a set of five coupled, non-linear differential 
equations describing the sheath. These are nondimensionalized to adjust all variables 
to the first order, and are rewritten below as equations (12)-(16), with 
nondimensionalized variables denoted by "x". Nondimensionalization can be 
accomplished as follows: The species number densities n„n„ are divided by their 
values at infinity to produce output from the anode surface to unity at the ambipolar 
boundary. The current densities j„j, are divided by the total current, allowing the 
output to show the "mirror behavior" of the two currents. The electric field is divided 
by the initial anode value to give output starting from unity at the surface and 
decreasing to the final core field value. The variable "y" is divided by the 
characteristic length’ of the plasma, "a", producing y. These corrections allow all 
output to vary in the range from zero to one, as a function of distance from the 
anode. 



’The characteristic length is defined so as to cancel the multiplying factor in the 
electric field equation, (14), (a = 1.107 x lO"). This allows a physical interpretation of 
the ion/electron number densities, as well as the decay rate of the electric field. 
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Attempts to solve this equation set using the computer code discussed below shows 
the set to be extremely sensitive to initial conditions. The computer code solver uses 
a "marching" scheme from the anode to the undisturbed plasma. The initial 
conditions are chosen to produce the electric field potential drop observed in actual 
thrusters. First and second space derivatives of the electric field are used as 
diagnostic checks to ensure reasonable output values and indicate instability of the 
integration process. Figure 4 shows the required resulting curves for the electric field 
and its first and second derivatives. 
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1.5 





Figure 4 - Electric Field and First Two Space Derivatives Used as Diagnostic Checks 
for Integrator Output. (Plotted for a Generic Exponential Function). 



19 



Ecker characterizes the plasma at the anode as a double sheath, with the inner 
section called the "inertia sheath", and an outer section called the "energy loss 
section". The inner section shows a potential rise of the order of one volt, with the 
outer section showing the exponential potential drop shown in Figure 4. While this 
double sheath may in fact describe the actual sheath region, the formulation above 
only models the potential drop portion of the sheath, and does not attempt to 
produce the potential rise of the irmer sheath. In addition, Ecker’s current 
constrictions are of a "macroscopic" nature, whereas those of Reference 14 and this 
work are "microscopic" [Ref. 33]. 

Data for a 6,000°K Nitrogen plasma were used to test the equation set [Ref. 
21]. Producing a proper solution required adjusting the initial conditions to force the 
curve shapes discussed above. Using Equations (2), (3), and (4), the mean free path, 
Debye length, and approximate sheath extent are calculated as X = 9 . 352xi0'^ m, 
Xp = 1.69 0x10*® m, and Xg = 2.352x10*® m (this assumes a drop voltage of 100 
volts). 
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3. Approximate Formulation 



Reference 21 explores the above equation set by taking advantage of the 
symmetry among the equations, and introducing two parameters, K' and K", shown 
below. 



K* = 



j, 

eD 

1 




(17) 



-K- 



j, . L 

e D eD 

I e 



(18) 



It can be shown that the resulting equations can be manipulated to yield a single, 
ordinary differential equation for the K’s in terms of the electric field. The resulting 
equation can be scrutinized for two distinct temperature regimes. Note that while 
the total current density, J, is constant in a steady, one-dimensional case, the K’s can 
vary and will in turn also depend on the degree of reactivity of the plasma (fi^), i.e.. 






= en 



(19) 



dy dy 

Because ion diffusion is much slower than electron diffusion, it can be shown that the 



K’s are related by 



K- = -K- + — (20) 

eD 

e 

As will be evident, at the electrode surface the K’s are equal to each other and at the 
undisturbed plasma, K' = 0. The total current density may be evaluated from 
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( 21 ) 



where v,. is the electron drift velocity beyond the ambipolar region which is strictly 
a function of E./N, (i.e., of the ratio of undisturbed electric field to the total number 
density). 

a. Effects of Temperature on Anode Constriction 

It is useful to investigate the overall effects of temperature. Since 
temperature will be considered constant, it comes in as a parameter in this 
formulation whereas charge density and electric field remain as variables. Intuitive 
arguments will be introduced which suggest that the electron and ion/neutral 
temperatures play a rather singular role in determining the intrinsic dimensionality 
of the problem, (i.e., there are cases when the geometry of the current lines is not 
necessarily impressed by the electrode geometry). Since the problem is described by 
moderate pressure, largely collisional sheaths, the ion and neutral temperatures are 
anticipated to remain reasonably equal. Depending on the gas, the electron 
temperature, on the other hand, can be elevated from the gas temperature at the 
anode where actual magnitudes depend on the local value of E/N. In order to get 
a perspective on the effects of temperature, we shall consider two extremes, namely, 
the case where the electron and ion temperature are the same (the equilibrium case) 
and the case where the electron temperature is substantially elevated from that of 
the ions/neutrals (the two-temperature case). 
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(1) Case I: T, = T, = (Equilibrium) 

The charge densities can be eliminated by combining equations 
(5)-(9), (17) and (18). The resulting equation can be shown to be 
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If the electric field decreases monotonically from the wall to the undisturbed plasma 
(i.e., from £„ - E.), then as y - «>, E E., E’ - 0, E” 0. 

So that in equation (22) above the "outer solution" becomes: 

K' = — (23) 

eD 

c 

Now this represents an acceptable solution from a physical point of view. Moreover, 
as y oo , 

n « D (K*)' = 0 (24) 

which is also acceptable for an equilibrium situation at the undisturbed plasma. 
Results [Ref. 21] are shown in Figure 5 for the case of nitrogen at 6000°K using an 
approximate electric field distribution. 
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Figure 5 - Electric Field and Species as a Function of y, Distance From Anode. An 
Approximation Using a Shaped Electric Field and Isothermal Plasma [Ref. 21). 
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(2) Case II: T, > > T, = T„ (Two -Temperature) 



In this case the same procedure as before yields the following 
equation where terms divided by T, have been dropped when compared to their 
counterparts divided by T,,. 

T 2 6 6 6 

(K*)' - _ E' = ^ + i EE" + _^E'E' - —E'" (25) 

E kT D kT eE e 

o e o 

Assuming the same monotonic decrease as before for the electric field from the wall 
to the plasma proper, as y -♦ oo, E -* E., E’ 0, E” -♦ 0. 

Then the outer solution becomes 



dK- ^ 2eEJ 
dy kT eD 

^ o e 



with h > 0 

e® 



(26) 



Or, K* -♦ (constant) y + constant, and keeps increasing with y. 

This is not the proper outer solution for the one-dimensional, equilibrium 
plasma that we seek because the net ionization rate continues to increase well inside 
the plasma proper where conditions should saturate, yielding a constant electric field. 
Therefore, as formulated. Case II is not amenable to a one-dimensional solution. 
References 14 and 21 show how this case can be analyzed under a multidimensional 
approach. These references also discuss a method for describing the electron 
temperature as a function of E/N, then how to couple a simplified energy relation 
which satisfactorily describes a two-temperature plasma. The necessary ingredient 
to make equation (26) approach zero beyond the decrease of E to E„ is to allow J 
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to fan out as indicated in Figure 6. Thus, in equation (26), the product "EJ" can 



bring down the charge production rate to arbitrarily low values. Alternatively, it is 



possible to explore techniques of bringing the electron temperature down to be in 



closer equilibrium with the ions and neutrals. Transpiration cooling is one such 



means. 



U(x,y) B J(x,y) 

c=C> 0 



Y 




Figure 6 - Two-Dimensional Model of Current Paths Showing Periodic Structure. 
Thermal Instabilities and Inhomogeneities Would Favor One Site Over Others and 
a Single Macroscopic Constriction May Then Be Produced. (Refs. 14. 34] 



26 



b. Similarity to Vacuum Arc Phenomena 

Instability phenomena observed in vacuum arcs [Ref. 35] are very 
similar to those observed in self-field thrusters [Ref. 12]. After the establishment of 
the current, the anode region operates in a vapor that issues from the electrodes. 
In vacuum arcs, Miller characterizes the anode region as operating in one of five 
distinct modes, ranging from a passive, low current mode to a high current, fully 
developed spot mode [Ref. 36]. Given the similarities mentioned above, vacuum arc 
anode research should be helpful in the understanding of MPD thruster transition 
to the anode spot mode. Existence diagrams after Miller [Ref. 36] are shown in 
Figure 7, which divide operating modes into regions as a function of anode current 
versus electrode geometry. Figure 7 shows the transition from glow to spot mode. 

Anode spot formation at high currents is clearly a factor in limiting anode 
lifetime. Various phenomena have been related to anode spotting. Hugel [Ref. 12] 
relates the transition to spotting mode to an increase in J^/m above a critical level. 
A separate factor connected with the spot mode is surface temperature of the anode. 
Rich, et.al., [Ref. 37] show that anode spotting is preceded by a luminous "footpoint" 
and followed by local melting prior to spot formation. Separately, Schuocker [Ref. 
38] finds a connection between spotting initiation and the factors of anode 
evaporation and magnetic constriction in vacuum arcs with high currents. 
Experimental investigations must be performed to see if the above-mentioned 
vacuum arc criteria apply to self-field thrusters. 
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Figure 7 - Anode Discharge Modes as a Function of Current and Gap Length. [Ref. 
36] 
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c. COMPUTER CODE 



Rather than using linear approximations to equations (12)-(16), the nonlinear 
set was used, with initial conditions adjusted in an attempt to produce observed 
electric fields from probe data. First and second space derivatives of the electric field 
were used as diagnostic checks to ensure computed output was reasonable. Initial 
conditions computed from the approximate formulae in Reference 21 were used. The 
equation set above presents a difficult problem for two reasons, nonlinearity and 
multiple time constants. The species number density equations, (12) and (13), both 
contain a nonlinear term, each with a time constant of its own. In addition, the 
electric field equation, (14), adds a possible third time constant. This constitutes a 
"stiff set of equations. Attempts were made to solve the set with the data discussed 
above, using Gear’s method of backward differentiation, in hopes that the variables 
would change slowly enough with each iteration to render a convergent iterative 
process. As described in Reference 39, if some reactions are slow and others fast 
among a set of coupled equations, the fast ones will control the stability of the 
method. This is addressed in the DGEAR program available from the International 
Mathematical & Statistical Library (IMSL). The latter software contains an Adams 
predictor-corrector method, as well as Gear’s method, which is well known for its 
success at solving stiff equation sets. The DGEAR software allows for a choice of 
functional or chord iteration methods, as well as a choice of Jacobian matrices. A 
more detailed discussion of this software can be found in Reference 39 and in the 
IMSL library. [Ref. 39] 
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D. COMPUTATIONAL RESULTS 



Numerous computer runs were completed using the initial conditions taken 
from Reference 21. In addition, data for the ionization coefficient v. Figure 8, was 
taken from References 40 and 41. 




Figure 8. Ionization Coefficient v as a Function of E/N for Nitrogen for Various 
Vibrational Temperatures (Refs. 40,41). 

Various combinations of initial conditions and ionization coefficient were used. As 
mentioned above, the electric field and its space derivatives were used as 
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diagnostic/reasonability checks on the output. Individual, as well as multiple 
computer runs were attempted to model the sheath region. Nonlinearities in the 
equation set are clearly seen in Figure 9. The ion number density does not reach that 
of the electrons, and the latter population growth rate continues to grow without 
bound. The shape of the electron population curve is very sensitive to its initial value. 
As shown in Figure 9, the latter population has too high a growth rate when 
compared to the ion population, and the latter does not "catch up". Increasing the 
initial value of fig flattens out this curve to a reasonable shape. Above an initial 
value of approximately 0.06, however, the plot of fig "dips" after a certain distance 
and then continues to increase as expected. This gives an approximate upper value 
for this initial value. To avoid instabilities like this, small "slices" were taken of the 
output after a small number of integration steps and multiple runs were used to form 
a "cut and paste" plot of the region. When a reasonable plot shape was produced, the 
value of ionization coefficient was varied in the "slices" to attempt to produce the 
required end values for electric field and species population. Both multiple and 
equilibrium values for the ionization coefficient were used. When the data showed 
signs of instability and failure to follow the required forms of Figure 4, a "slice" was 
made in the data stream, and the data points from this point used to start a new 
computer run. This approach was taken in the hope of avoiding singularities in the 
integration from anode surface to ambipolar region. In addition to the diagnostic 
checks shown in Figure 4, an additional data check is provided by the transition from 
the sheath to the ambipolar region. 
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Figure 9. Species Number Density plots for Individual Computer Run, Showing 
Divergent Tracks for Ion and Electron Populations, and Effect of Nonlinearity. 



As shown in Figure 5, the species number densities are equivalent in this region, as 
are their change rate. Thus, setting Equations (12) and (15) equal to each other and 
solving for n yields a value of 0.5 in the ambipolar region. As indicated in Figures 
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10-11, the output produces the desired plot slopes for electric field and species 
number density. However, the number density plots cross long before approaching 
the required value of 0.5. In addition, neither electric field nor species number 
density approaches an equilibrium value or shows sign of levelling off. Apparently, 
the multiple time constants and nonlinear portions of the number density equations 
combine to create a seemingly intractable system. Solutions for this system may be 
possible for specific, individual initial condition sets, but the problem does not appear 
amenable to this approach in general. A one-dimensional system such as this may be 
better described through the approach of boundary layer theory or nonlinear 
dynamics and chaos. Given the effort and difficulty involved in the latter, a one- 
dimensional approach such as that modelled above does not appear useful. A 
combination of one- and two-dimensional modelling would appear to be more useful, 
as discussed in Reference 14. A one-dimensional model may be useful, but only in 
an approximation approach, with a shaped electric field, such as that used in 
Reference 21. 
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Figure 10. Electric Field as a Function of y. Distance From Anode, Using 
Equations ( 12)-(16). 
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Figure 1 1. Species Number Density as a Function of y , Distance From Anode, Using 
Equations (12)-(16). 
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E. ANODE FALL VOLTAGE 



The arc discharge operates with an anode voltage drop (V„) which should be 
on the order of the ionization potential of the gas as it is in glow discharges. This 
voltage drop is largely non-ohnric because the anode region is typically very thin and 
space-charge controlled. Since a significant portion of the anode heating comes via 
the oncoming electrons accelerating through the fall voltage, it is of interest to 
minimize this voltage, which can range from less than ten to over 40 volts in argon. 
The question arises as to what governs this anode drop and why is it such a 
noticeable function of current? [Ref. 42] 

In collisional sheaths, the anode voltage drop is largely governed by a positive 
ion generation region which forms in front of the anode. The production of ions 
reduces the space charge density and thereby permits operation at lower voltages 
than otherwise possible. Ions are most often created by electrons within their last 
few mean-free-paths before entering the anode; these ions slowly drift away from the 
anode, thereby effectively neutralizing the space charge, at a speed proportional to 
the ratio of their mobility to that of the electrons. At the anode, the electric field 
is a maximum and the electron mean energy displays a corresponding high. 

At moderate pressures, the sheath is very thin and breakdown can be visualized 
as occurring between the undisturbed plasma (which acts as a rich source of 
electrons) and the anode surface. This arrangement has the attributes that represent 
"thermionic arc breakdown" [Ref. 43], a source of electrons which is independent of 
the breakdown itself and relatively small spacings. For gases that allow cumulative 
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ionization (with some help from the tail of the Maxwellian distribution of electron 
energies), what results is a breakdown voltage appreciably below the ionization 
potential. This then could be an explanation for the low voltage breakdown 
observations [Ref. 42]. Clearly, gases with low ionization potentials and lots of 
atomic electron energy levels are preferred (such as cesium and barium) but low- 
voltage breakdown has been observed with most gases. 

The increase of the anode fall voltage above the ionization potential has been 
related to the electron Hall parameter, since a reduction of this parameter decreases 
that voltage and corresponding losses. Control of the local magnetic field through 
the use of and array of permanent magnets as well as implementation of 
transpiration cooling (which increases the electron collision frequency) have both 
yielded some encouraging results. Because the anode fall also scales up with JVA> 
it is conceivable that current inhomogeneities and plasma instabilities which are 
reflected in this parameter are in the picture as well. [Ref. 44] 

In summary, any possible reduction of the anode fall voltage will hinge on a 
thorough understanding of the anode region, with its associated sheath and ambipolar 
regions, where electron temperature effects, ionization effects, and magnetic field 
effects play a pivotal role. If transpiration cooling is present, then additional 
phenomena of fluid-dynamic nature may come into play. Experimental observations 
with atmospheric discharges indicate the possible presence of convective effects at 
the anode. [Ref. 45] 
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IV. TRANSPIRATION COOLING 



Transpiration cooling of the anode has often been promoted as an attractive 
means of recovering a large portion of the power deposited there. Additionally, the 
onset of melting may be minimized or even avoided by active anode cooling. Rich, 
et.al., related high anode temperatures to anode spotting [Ref. 37]. Similarly, Park 
and Choi showed that low thermal diffusion leads to erosion and, consequently, 
anode damage [Ref. 46]. Active anode cooling via transpiration is one means of 
ensuring high thermal diffusion and extending anode lifetime. Early work by 
Schoeck, et. al., [Ref. 47] showed that up to 80% of the energy deposited in the 
anode is recoverable via transpiration cooling. While this study used non-convective, 
high-intensity argon arcs, it is reasonable to assume that this effect would apply in 
MPD arcs using other propellants. Although this cooling method has not been 
studied for incorporation in plasma thrusters for some time, it has been recently 
considered as a means of cooling the fuselage of the National Aerospace Plane 
(NASP). Plasma thruster designs could undoubtedly gain from this database, and due 
consideration should be given to this cooling approach for the anode. 

For a given mass transfer flow rate, the heat flux reduction to a surface is 
inversely proportional to the molecular weight of the injected gas. Use of the 
propellant as coolant as well as fuel would eliminate the need for additional tankage 
and pumps, simplifying the design considerably. Lithium has been considered to be 
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the propellant of choice, primarily because of its low molecular weight, its favorable 

ionization potential, and its low-volume tankage properties. It has a relatively low 

first ionization potential of 5.4 eV and a high second ionization potential of 75.6 eV. 

This single ionization potential range of over 70 eV compares to approximately 20 

eV for Cesium and 27 eV for Potassium [Ref. 48]. This provides a broad temperature 

range within which only single ionization will occur. Large temperatures must be 

reached within the gas before double ionization occurs. As the gas temperature is 

increased several thousand degrees Kelvin, it undergoes ionization and disassociation. 

Thermal energy deposited can be recovered through nozzle expansion at the exit. 

However, residence time of the gas is not long enough to ensure recombination. 

Thus, the energy invested in ionization and disassociation will be lost. [Ref. 49] 

Lithium has been shown to produce specific impulse figures in excess of 7,000 

seconds at 70% efficiency in steady-state thrusters, [Ref. 50] whereas all other 

propellants have been limited to less than 3,000 seconds specific impulse at less than 

40% efficiency [Ref. 6]. Subramaniam has concluded that: 

...regenerative cooling of anodes (at the specific impulse values in the MPD 
regime) is possible only with hydrogen or with alkali-metal propellants, notably 
lithium. In the latter case, the ideal anode operating mode would be 
evaporation and ionization of the propellant on the porous or wetted anode 
surface, resulting in increased ion current fraction, reduced anode voltage fall 
and utilization of part of the anode loss energy [Ref. 51]. 

Liquid coolants, as well as reducing storage requirements, offer the advantage of 

providing latent heat of vaporization for energy disposal. However, design problems 

can occur if the liquid is allowed to vaporize within the porous structure. Problems 
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arise due to the abrupt increase in pressure gradient as the coolant vaporizes. Since 
coolant flow generally have three-dimensional characteristics, the flow will be 
diverted around the vapor bubbles and hot spots often develop. The technical 
practicality of using molten lithium to cool a porous tungsten anode would seem to 
be beyond current technology. On the other hand, the products of decomposition of 
hydrazine (gaseous hydrogen and ammonia) have proven to be efficient and practical 
coolants [Ref. 52]. 

Given the performance figures above, using an auxiliary coolant gas even with 
high molecular weight (e.g., NH 3 , N^, CH 4 , etc.) which could serve as a propellant 
once released from a porous hot tungsten anode surface would seem to more 
practical, vice dealing with molten lithium. Experimental studies would be needed 
to compare the approaches. Kuriki and Suzuki performed experiments with a quasi- 
steady MPD thruster to study the effect of anode gas injection (Argon). At high 
currents of up to 10 kA, increases in thrust, specific impulse, and flow discharge 
stability were observed [Ref. 53]. 

There is some question as to the likelihood of current constriction resulting 
from anode gas injection. In such a case, swirling or circulating the propellant gas 
would help to move any footpoints that developed around the anode surface and 
prevent them from becoming fully-developed spots. Additionally, an applied 
magnetic field could serve to circulate the footpoints as well. The unique advantage 
of transpiration cooling hinges on providing effective anode cooling while supplying 
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hot propellant, but the real benefit will depend on how small the amount of coolant 
required really will be. 

Transpiration cooling has proven to be as desirable as it is challenging. It is 
complicated to implement, with associated reliability problems and difficulty of 
analytical predictions. While the production of thicker boundary layers is largely 
ineffective against the electron flux heating, the cooling itself is most efficient and a 
substantial fraction of the energy transferred to the anode is recoverable. The 
arguments of Chapter Three indicate that a reduction of the electron temperature 
in the anode would have the desirable effect of reducing the initial current spotting 
which can be conjectured to be the path that leads to anode arc spots. This electron 
temperature reduction can be done most effectively by polyatomic gases (which have 
a high 5-loss factor) emanating from the anode surface [Ref. 54]. 

The arguments relating to transpiration cooling might be summarized as 
follows; 

Favorable Outcomes 

• No separate cooling mechanism for anode required, 

• Adds "hot" propellant to exhaust "recovering" most of the electrical power loss to 
the anode, 

• Quenches T^ thus likely to postpone anode spotting and reduce the heating 
associated with the electron thermal energy (5kT,/2e), 

• Reduces bulk convective heating, 

• Reduces the local electron Hall parameter by increasing the collision frequency. 
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Favorable Outcomes (cont’d.) 

• Allows for some radiation cooling from the hot tungsten surface (about 120 
watts/cm‘ at 2800°K [Ref. 51]), 

• Hydrogen/ammonia gases flowing through hot porous (sintered) tungsten 
represent a compatible, proven technology. 

Unfavorable Outcomes 

• May decrease the electrical conductivity in the anode region, 

• May destabilize the ionization processes in the sheath and bring about significant 
fluctuations in the current, 

• Disrupts "cathode jet" in front of the anode with unpredictable consequences, 

• Introduces propellant which may not be hot enough, not ionized enough, or not 
in the proper place for J xB acceleration, 

• Transpiration cooling through a porous (tungsten) anode is a difficult design 
problem. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



Plasma thrusters offer distinct advantages in terms of payload delivered for 
interplanetary missions, as well as for orbital transfer. A recent comparison 
completed by Choueiri, Kelly, and Jahn shows a mass savings of 65 tons for an 
orbital transfer from low Earth orbit to geosynchronous Earth orbit using a 
quasisteady MPD thruster as opposed to an advanced chemical thruster.^ This 
superior performance comes at the expense of low thrust-to-weight ratio and long 
transit time. However, given the large cargo/logistic requirements of a manned 
interplanetary mission, delivery of payload must be maximized. Thus, further work 
to characterize more fully thruster behavior and anode contributions in particular are 
certainly warranted. [Ref. 55] 

The "cut-and-paste" method used to generate Figures 10 and 11 is not of 
practical use as a modelling method, due to the large effort involved. It did produce 
the expected electric field and species number and current density plots near the 
anode, but failed to produce the entire sheath out to the ambipolar region. The 
nonlinearity of the equation set led to a quickly deteriorating solution. A more 
practical approach using nonlinear dynamics and/or chaos must be developed to 
model the sheath numerically. 



*This assumes a specific impulse of 2,000 seconds, 600 kW of input power, and 
a 270-day transit time. 
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Depending on the propellant mass fraction used for cooling, the transpiration 
scheme discussed above presents some rather unique advantages. A hot anode 
which uses only a small amount of propellant for cooling need not be penalized for 
any lost thrust. If in addition, we increase anode lifetime by delaying the formation 
of anode arc spots, then the scheme is all the more desirable. A decrease of the 
electron temperature in the vicinity of the anode may bring about a more 
homogeneous flow of current and a reduction in the heating effect associated with 
the high electron kinetic energy. Recovery of the heat deposited at the anode would 
be most important if the propellant fraction is high. In such case, nozzle expansion 
of the hot-propellant/coolant-gas might be implemented. 

Means of limiting anode losses through decreasing anode fall voltages were 
discussed, including the control of the local Hall parameter and the implementation 
of thermionic arc breakdown. The electrical conductivity (of a nonreacting plasma) 
could possibly decrease as a result of transpiration cooling and this might increase 
the anode fall voltage. 

Additional work needs to be done in the following areas: 

♦ Investigate effectiveness of nonlinear dynamics and chaos in solving sheath 
equation set, 

• Incorporate adequate one- or two-dimensional sheath modelling in quasisteady 
MPD numerical codes. 
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• Investigate the role that fluid dynamic effects play in MPD thruster anode 
discharges, 

• Investigate the effect of transpiration cooling on current and plasma stability, as 
well as on thruster performance and lifetime, 

• Determine effectiveness of transpiration cooling’s increase of the collision 
frequency parameter, 

• Compare performance of gaseous propellant/coolants versus hybrid designs with 
lithium propellant/gaseous coolant, 

• Determine if required percentage of propellant gas as coolant is practical (e.g., less 
than 10%), 

• Investigate effect of surface imperfections as focal points for current constrictions 
and as precursors to anode spotting. 
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APPENDIX A 



The following software includes the calling program, SHEATH, its two 
subroutines, FCNJ and YDOT, and the DGEAR integrator. The latter is quite 
extensive in length and includes ten subroutines, including the following: DGRST, 
DGRCS, DGRPS, DGRIN, LUDATF, LUELMF, LEQTIB, UERTST, UGETIO, and USPKD . A 
detailed discussion may be found in IMSL literature or Reference 39 . 

Program Sheath 000010 

C 000020 

C Calling program for DGEAR integrator. Initial conditions are 000030 

C input via READ statements and keyboard entry. Output is to data 000040 

C files via the DGRST subroutine. Diagnostic check of output via 000050 

C Figure 4 printed to data file from DGRST subroutine. Consult 000060 

C DGEAR comments for variable descriptions not listed below. 000061 

C--- 000062 

REAL E,K,EPS,TI,EF,EFINF,DI,DE,NUINF,C1,K1,A 000070 

INTEGER N, METH, MITER, INDEX, IWK(l) , lER, STEP 000080 

REAL*8 X,H, Y(5) ,XEND,TOL,WK 000090 

EXTERNAL YDOT, FCNJ, DGEAR 000100 

COMMON/CONST/E, K, EPS, TI,EF,EFINF, NINE, DI,DE, VINE, Cl, Kl, A 0 00110 

C 000120 

C Constants 000121 

C Cl and Kl are constants describing the ionization coefficient. 000122 

C They are taken from the data plotted in Figure 8. The 000123 

C coefficient is equal to the nondimensionalized electric potential 000124 

C raised to the Kl power and then multiplied by Cl. 000125 

C In this way, the ionization coefficient is allowed to vary in 000126 

C proportion to the strength of the electric potential. 000127 

C 000128 

WRITE (*,*)' Input value for Cl (format 6E3):' 000129 

READ(*,*)C1 000130 

WRITE (*,*) Cl 000131 

WRITE (*, ^) ' Input value for Kl (format 6E3):' 000132 

READ(*,*)K1 000133 

WRITE (*,*) Kl 000134 

C Initial conditions for species number density, electric potential 000135 

C and species current density are now input (ni , ne , E, je , j i) . 000136 

C - 000137 

WRITE (*,*)' Input values for y(l) through y(5) (format 5(6E3)):' 000138 

READ(*,*)y(l) ,y (2) ,y (3) ,y (4) ,y (5) 000139 

WRITE (*, *)y(l) ,y(2) ,y(3) ,y(4) ,y(5) 000140 

C Following constants are for plasma described in Reference 21 000141 

C (6,000 K, Init E=20,000 V/m, Final E=l,200 V/m) 000142 

E=1.6E-19 000143 

K=1.38E-23 000150 

EPS=8 .854E-12 000160 

TI=6E3 000170 

EF=2E5 000180 

EFINF=1 .2E4 000190 

DI=1 .724E-4 000200 

DE=1 .724E-1 000210 

VIO = 2.E6 000215 

VINE = 4.93E-7 000220 

C 000230 
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C A is plasma characteristic length which shows potential drop. 

C A= { (EPS*EF) / (E*NINF) ) = 1.107E-6 

A = 1.107E-5 

X = 0.01 

XEND = 10. 

H = le-6 
TOL = IE -6 
METH = 2 
MITER = 1 
INDEX = 1 
N=5 

IWK(l) = 5 
WK = 18000. 
lER = 0 

OPEN (UNIT=8,FILE=' SHEATH. DAT' , STATUS =' UNKNOWN' ) 

CALL DGEAR2 (N, YDOT, FCNJ, X, H, Y, XEND, TOL, METH, MITER, INDEX, IWK,WK, 
+IER,STEP) 

DO 3 1=0, N 
DO 2 J=0,100 

WRITE{*,*) J,Y(I) 

WRITE(8,1) J,Y(I) 

1 FORMAT{T2,F5.1,5 (5X,D9.2) ) 

2 CONTINUE 

3 CONTINUE 

WRITE {*,*) 'Total Steps = ', STEP, ' Final Step Size = ',H, 

+' Error Code = ' , lER 
CLOSE (UNIT=8) 

STOP 

END 

Q’k’kieie*************-k**’*r-k-k** ********** 

C DUMMY SUBROUTINE FCNJ 

Q** ******************************** * 

SUBROUTINE FCNJ {N, X, Y, PD) 

INTEGER N 

REAL Y(N) , PD (N,N) ,X 

RETURN 

END 

Q** ******************************* ** 

C SUBROUTINE YDOT 

Q****** ******************* ********** 

SUBROUTINE YDOT (N, X , Y , YPRIME , epr ime , epr ime2 ) 

REAL* 8 X, Y{5) , YPRIME (5) ,NUI,eprime,eprime2 
REAL E , K, EPS , TI , EF , EFINF , NINF , DI , DE , VINE , Cl , K1 , A, B1 , B2 , B3 , B4 
COMMON/CONST/E , K, EPS , TI , EF , EFINF , NINF , DI , DE , VIO , VINE , Cl , K1 , A 
VI = Cl * (Y(3)**K1) 

VIT = VI / VIO 

Following constants are the bracketed values in Equations 12-16. 
A is left as a variable. 



000240 

000250 

000270 

000280 

000290 

000300 

000305 

000310 

000320 

000330 

000340 

000350 

000360 

000370 

000380 

000390 

000400 

000410 

000420 

000430 

000440 

000450 

000460 

000470 

000480 

000490 

000500 

000510 



1 

2 

3 

4 

5 



B1 = ( (E*EPS) / (K*TI) ) * A 

Bl= 3.86E5 * A 

B2 = ( (E*EFINF) / (K*TI) ) * A 

B2 = 2.32E4 * A 

B3 = ( (E*NINF) / (EF*EPS) ) * A 

B3 = 9 . 04E5 * A 

B4 = ( {VINF*K*TI) / (E*DE*EFINF) ) * A 

B4 = 2.62E-21 * A 

Alpha = 2 -body recombination coefficient (fm. Laser Kinetics 
Handbook (AFWL-TR- 74 -216 , 1974)) (cm3/sec) 

Alpha = 9.e-8 
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non o o n n no 



C-- 

C FIVE FIRST ORDER EQUATIONS - Equations 12-16 



Ni 

YPRIME(l) = 


(B 


★ 


Y(l) * Y(3) ) - 


Y(5) 




Ne 

YPRIME(2) = 


- (B 


★ 


Y(2) * Y(3)) + 


Y(4) 




E 

YPRIMEO) = 


B3 


★ 


(Y(l) - Y(2)) 






je 

YPRIME(4) = 


-B4 


★ 


Y(2) * (VIT - 


(ALPHA * 


Y(l) ) ) 


ji 

YPRIME(5) = 


B4 


★ 


Y(2) * (VIT - 


(ALPHA * 


Y(l) ) ) 



Diagnostic Check of first, second derivatives 

eprime = y(l) - y(2) 

eprime2 = yprime(l) - yprime(2) 

C 

RETURN 

END 
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IMSL ROUTINE NAME 
-modified to return 
COMPUTER 
LATEST REVISION 
PURPOSE 

USAGE 

ARGUMENTS N 
FCN 



DGEAR DGEAOOlO 

DGEA0020 

# of steps via variable "step” in subroutine call + 

DGEA0040 
DGEA0050 
DGEA0060 
DGEA0070 
DGEA0080 

- DIFFERENTIAL EQUATION SOLVER - VARIABLE ORDER DGEA009 0 



of steps via 

- IBM/DOUBLE 

- NOVEMBER 1, 1984 



FCNJ 



DGEAOlOO 
DGEAOllO 
DGEA0120 
DGEA0130 
DGEA0140 
DGEA0150 
DGEA0160 
DGEA0170 
DGEA0180 
DGEA0190 
DGEA0200 
DGEA0210 
DGEA0220 
DGEA0230 
DGEA0240 
DGEA0250 
DGEA0260 
DGEA0270 
DGEA0280 
DGEA0290 
DGEA0300 
DGEA0310 
DGEA0320 
DGEA0330 
DGEA0340 
DGEA0350 
DGEA0360 
DGEA0370 
DGEA0380 
DGEA0390 
DGEA0400 
DGEA0410 
DGEA0420 
DGEA0430 
DGEA0440 
DGEA0450 
DGEA0460 
DGEA0470 
DGEA0480 
DGEA0490 
DGEA0500 
DGEA0510 
DGEA0520 
DGEA0530 
DGEA0540 

FCNJ MUST EVALUATE PD IN BAND STORAGE MODE. DGEA0550 
THAT IS, PD (N* (J-I+NLC) +1) IS THE PARTIAL DGEA0560 
DERIVATIVE OF YPRIME(I) WITH RESPECT TO DGEA0570 
Y(J) . NLC IS THE NUMBER OF LOWER DGEA0580 

CODIAGONALS FOR THE BAND MATRIX. DGEA0590 

FCNJ MUST APPEAR IN AN EXTERNAL STATEMENT INDGEA0600 
THE CALLING PROGRAM AND N, X, Y (1) , . . . , Y (N) DGEA0610 



ADAMS PREDICTOR CORRECTOR METHOD OR 
GEARS METHOD 

- CALL DGEAR (N, FCN, FCNJ, X, H, Y, XEND , TOL, METH , 

MITER, INDEX, IWK,WK, lER) 

- INPUT NUMBER OF FIRST-ORDER DIFFERENTIAL 

EQUATIONS . 

- NAME OF SUBROUTINE FOR EVALUATING FUNCTIONS. 

(INPUT) 

THE SUBROUTINE ITSELF MUST ALSO BE PROVIDED 
BY THE USER AND IT SHOULD BE OF THE 
FOLLOWING FORM 

SUBROUTINE FCN (N,X, Y, YPRIME) 

REAL X,Y(N) , YPRIME (N) 



FCN SHOULD EVALUATE YPRIME (1) , . . . , YPRIME (N) 
GIVEN N,X, AND Y(l) , . . . ,Y(N) . YPRIME (I) 

IS THE FIRST DERIVATIVE OF Y(I) WITH 
RESPECT TO X. 

FCN MUST APPEAR IN AN EXTERNAL STATEMENT IN 
THE CALLING PROGRAM AND N, X, Y (1 ),..., Y (N) 
MUST NOT BE ALTERED BY FCN. 

- NAME OF THE SUBROUTINE FOR COMPUTING THE 

JACOBIAN MATRIX OF PARTIAL DERIVATIVES. 

(INPUT) 

THE SUBROUTINE ITSELF MUST ALSO BE PROVIDED 
BY THE USER. 

IF MITER=1 IT SHOULD BE OF THE FOLLOWING 
FORM 

SUBROUTINE FCNJ (N,X,Y,PD) 

REAL X, Y(N) ,PD (N,N) 



FCNJ MUST EVALUATE PD(I,J) , THE PARTIAL 
DERIVATIVE OF YPRIME (I) WITH RESPECT TO 
Y(J) , FOR 1=1, N AND J=1,N. 

IF MITER= -1 IT SHOULD BE OF THE FOLLOWING 
FORM 

SUBROUTINE FCNJ (N,X,Y,PD) 

REAL X,Y(N) ,PD(1) 
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H 



XEND 



TOL 



METH 



MITER 



MUST NOT BE ALTERED BY FCNJ . DGEA062 0 

FCNJ IS USED ONLY IF MITER IS EQUAL TO DGEA0630 

1 OR -1. OTHERWISE A DUMMY ROUTINE CAN DGEA0640 
BE SUBSTITUTED. SEE REMARK 1. DGEA0650 

- INDEPENDENT VARIABLE. (INPUT AND OUTPUT) DGEA0660 

ON INPUT, X SUPPLIES THE INITIAL VALUE DGEA0670 

AND IS USED ONLY ON THE FIRST CALL. DGEA0680 

ON OUTPUT, X IS REPLACED WITH THE CURRENT DGEA0690 
VALUE OF THE INDEPENDENT VARIABLE AT WHICHDGEA0700 
INTEGRATION HAS BEEN COMPLETED. DGEA0710 

- INPUT /OUTPUT. DGEA0720 

ON INPUT, H CONTAINS THE NEXT STEP SIZE IN DGEA0730 
X. H IS USED ONLY ON THE FIRST CALL. DGEA074 0 

ON OUTPUT, H CONTAINS THE STEP SIZE USED DGEA0750 
LAST, WHETHER SUCCESSFULLY OR NOT. DGEA0760 

- DEPENDENT VARIABLES, VECTOR OF LENGTH N. DGEA0770 

(INPUT AND OUTPUT) DGEA0780 

ON INPUT, Y(l) , . . . , Y(N) SUPPLY INITIAL DGEA0790 

VALUES. DGEA0800 

ON OUTPUT, Y(l) , . . . , Y(N) ARE REPLACED WITH DGEA0810 
A COMPUTED VALUE AT XEND. DGEA0820 

- INPUT VALUE OF X AT WHICH SOLUTION IS DESIRED DGEA0830 



NEXT. INTEGRATION WILL NORMALLY GO 



DGEA0840 



BEYOND XEND AND THE ROUTINE WILL INTERPOLATEDGEA0850 
TO X = XEND. DGEA0860 

NOTE THAT (X-XEND) *H MUST BE LESS THAN DGEA0870 

ZERO (X AND H AS SPECIFIED ON INPUT) . DGEA0880 

- INPUT RELATIVE ERROR BOUND. TOL MUST BE DGEA0890 

GREATER THAN ZERO. TOL IS USED ONLY ON THE DGEA0900 
FIRST CALL UNLESS INDEX IS EQUAL TO -1. DGEA0910 

TOL SHOULD BE AT LEAST AN ORDER OF DGEA0920 

MAGNITUDE LARGER THAN THE UNIT ROUNDOFF DGEA0930 

BUT GENERALLY NOT LARGER THAN .001. DGEA0940 

SINGLE STEP ERROR ESTIMATES DIVIDED BY DGEA0950 

YMAX(I) WILL BE KEPT LESS THAN TOL IN DGEA0960 

ROOT -MEAN -SQUARE NORM (EUCLIDEAN NORM DGEA0970 

DIVIDED BY SQRT(N)) . THE VECTOR YMAX OF DGEA0980 
WEIGHTS IS COMPUTED INTERNALLY AND STORED DGEA0990 
IN WORK VECTOR WK. INITIALLY YMAX (I) IS DGEAIOOO 

THE ABSOLUTE VALUE OF Y(I) , WITH A DEFAULT DGEAIOIO 
VALUE OF ONE IF Y(I) IS EQUAL TO ZERO. DGEA1020 

THEREAFTER, YMAX (I) IS THE LARGEST VALUE DGEA1030 
OF THE ABSOLUTE VALUE OF Y(I) SEEN SO FAR, DGEA1040 
OR THE INITIAL VALUE OF YMAX (I) IF THAT IS DGEA1050 
LARGER. DGEA1060 

- INPUT BASIC METHOD INDICATOR. DGEA1070 

USED ONLY ON THE FIRST CALL UNLESS INDEX IS DGEA1080 
EQUAL TO -1. DGEA1090 

METH = 1, IMPLIES THAT THE ADAMS METHOD IS DGEAllOO 
TO BE USED. DGEAlllO 

METH = 2, IMPLIES THAT THE STIFF METHODS OF DGEA1120 
GEAR, OR THE BACKWARD DIFFERENTIATION DGEA1130 

FORMULAE ARE TO BE USED. DGEA114 0 

- INPUT ITERATION METHOD INDICATOR. DGEA1150 

MITER = 0, IMPLIES THAT FUNCTIONAL DGEA1160 

ITERATION IS USED. NO PARTIAL DGEA1170 

DERIVATIVES ARE NEEDED. A DUMMY FCNJ DGEA1180 

CAN BE USED. DGEA1190 

MITER = 1, IMPLIES THAT THE CHORD METHOD DGEA1200 
IS USED WITH AN ANALYTIC JACOBIAN. FOR DGEA1210 
THIS METHOD, THE USER SUPPLIES DGEA1220 

SUBROUTINE FCNJ. DGEA1230 
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INDEX 



IWK 

WK 



MITER = 2, IMPLIES THAT THE CHORD METHOD DGEA124 0 

IS USED WITH THE JACOBIAN CALCULATED DGEA1250 

INTERNALLY BY FINITE DIFFERENCES. DGEA1260 

A DUMMY FCNJ CAN BE USED. DGEA1270 

MITER = 3, IMPLIES THAT THE CHORD METHOD DGEA12 80 

IS USED WITH THE JACOBIAN REPLACED BY DGEA1290 

A DIAGONAL APPROXIMATION BASED ON A DGEA1300 

DIRECTIONAL DERIVATIVE. DGEA1310 

A DUMMY FCNJ CAN BE USED. DGEA132 0 

MITER = -1 OR -2, IMPLIES USE THE SAME DGEA1330 

METHOD AS FOR MITER= 1 OR 2 , RESPECTIVELY , DGEA13 40 
BUT USING A BANDED JACOBIAN MATRIX. IN DGEA135 0 
THESE TWO CASES BANDWIDTH INFORMATION DGEA1360 

MUST BE PASSED TO DGEAR THROUGH THE DGEA1370 

COMMON BLOCK DGEA1380 

COMMON /DBAND/ NLC,NUC DGEA1390 

WHERE NLC=NUMBER OF LOWER CODIAGONALS DGEA14 00 

NUC=NUMBER OF UPPER CODIAGONALS DGEA1410 

INPUT AND OUTPUT PARAMETER USED TO INDICATE DGEA1420 
THE TYPE OF CALL TO THE SUBROUTINE. ON DGEA1430 

OUTPUT INDEX IS RESET TO 0 IF INTEGRATION DGEA1440 
WAS SUCCESSFUL. OTHERWISE, THE VALUE OF DGEA1450 

INDEX IS UNCHANGED. DGEA1460 

ON INPUT, INDEX = 1, IMPLIES THAT THIS IS THE DGEA1470 
FIRST CALL FOR THIS PROBLEM. DGEA1480 

ON INPUT, INDEX = 0, IMPLIES THAT THIS IS NOT DGEA1490 
THE FIRST CALL FOR THIS PROBLEM. DGEA1500 

ON INPUT, INDEX = -1, IMPLIES THAT THIS IS NOTDGEA1510 
THE FIRST CALL FOR THIS PROBLEM, AND THE DGEA1520 
USER HAS RESET TOL . DGEA1530 

ON INPUT, INDEX = 2, IMPLIES THAT THIS IS NOT DGEA1540 
THE FIRST CALL FOR THIS PROBLEM. INTEGRATIONDGEA1550 
IS TO CONTINUE AND XEND IS TO BE HIT EXACTLYDGEA1560 
(NO INTERPOLATION IS DONE) . THIS VALUE OF DGEA1570 
INDEX ASSUMES THAT XEND IS BEYOND THE DGEA1580 

CURRENT VALUE OF X. DGEA159 0 

ON INPUT, INDEX = 3, IMPLIES THAT THIS IS NOT DGEA1600 
THE FIRST CALL FOR THIS PROBLEM. INTEGRATIONDGEA1610 
IS TO CONTINUE AND CONTROL IS TO BE RETURNEDDGEAl 6 2 0 
TO THE CALLING PROGRAM AFTER ONE STEP. XEND DGEA1630 



IS IGNORED. 



DGEA1640 



INTEGER WORK VECTOR OF LENGTH N. USED ONLY IF DGEA1650 
MITER = 1 OR 2 DGEA1660 

REAL WORK VECTOR OF LENGTH 4^N+NMETH+NMITER . DGEA1670 

THE VALUE OF NMETH DEPENDS ON THE VALUE OF DGEA1680 
METH. DGEA1690 

IF METH IS EQUAL TO 1, DGEA1700 

NMETH IS EQUAL TO N*13 . DGEA1710 

IF METH IS EQUAL TO 2, DGEA1720 

NMETH IS EQUAL TO N*6. DGEA1730 

THE VALUE OF NMITER DEPENDS ON THE VALUE OF DGEA1740 
MITER. DGEA1750 

IF MITER IS EQUAL TO 1 OR 2 , DGEA1760 

NMITER IS EQUAL TO N* (N+1) DGEA1770 

IF MITER IS EQUAL TO -1 OR -2, DGEA1780 

NMITER IS EQUAL TO (2 *NLC+NUC+3 ) *N DGEA1790 

WHERE NLC=NUMBER OF LOWER CODIAGONALS DGEA1800 
NUC=NUMBER OF UPPER CODIAGONALS DGEA1810 
IF MITER IS EQUAL TO 3, DGEA1820 

NMITER IS EQUAL TO N. DGEA1830 

IF MITER IS EQUAL TO 0, DGEA1840 

NMITER IS EQUAL TO 1. DGEA1850 
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lER 



WK MUST REMAIN UNCHANGED BETWEEN SUCCESSIVE 
CALLS DURING INTEGRATION. 

ERROR PARAMETER. (OUTPUT) 

WARNING ERROR 

lER = 33, IMPLIES THAT X+H WILL EQUAL X ON 
THE NEXT STEP. THIS CONDITION DOES NOT 
FORCE THE ROUTINE TO HALT. HOWEVER, IT 
DOES INDICATE ONE OF TWO CONDITIONS. 

THE USER MIGHT BE REQUIRING TOO MUCH 
ACCURACY VIA THE INPUT PARAMETER TOL . 

IN THIS CASE THE USER SHOULD CONSIDER 
INCREASING THE VALUE OF TOL. THE OTHER 
CONDITION WHICH MIGHT GIVE RISE TO THIS 
ERROR MESSAGE IS THAT THE SYSTEM OF 
DIFFERENTIAL EQUATIONS BEING SOLVED 
IS STIFF (EITHER IN GENERAL OR OVER 
THE SUBINTERVAL OF THE PROBLEM BEING 
SOLVED AT THE TIME OF THE ERROR) . IN 
THIS CASE THE USER SHOULD CONSIDER 
USING A NONZERO VALUE FOR THE INPUT 
PARAMETER MITER. 

WARNING WITH FIX ERROR 

lER = 66, IMPLIES THAT THE ERROR TEST 
FAILED. H WAS REDUCED BY .1 ONE OR MORE 
TIMES AND THE STEP WAS TRIED AGAIN 
SUCCESSFULLY. 

lER = 67, IMPLIES THAT CORRECTOR 
CONVERGENCE COULD NOT BE ACHIEVED. 



DGEA1860 
DGEA1870 
DGEA1880 
DGEA1890 
iJGEAlSOO 
DGEA1910 
DGEA1920 
DGEA1930 
DGEA1940 
DGEA1950 
DGEA1960 
DGEA1970 
DGEA1980 
DGEA1990 
DGEA2000 
DGEA2010 
DGEA2020 
DGEA2030 
DGEA2040 
DGEA2050 
DGEA2060 
DGEA2070 
DGEA2080 
DGEA2090 
DGEA2100 
DGEA2110 
DGEA2120 
DGEA2130 

H WAS REDUCED BY .1 ONE OR MORE TIMES AND DGEA214 0 
THE STEP WAS TRIED AGAIN SUCCESSFULLY. DGEA2150 
TERMINAL ERROR DGEA2160 

lER = 132, IMPLIES THE INTEGRATION WAS DGEA2170 

HALTED AFTER FAILING TO PASS THE ERROR DGEA2180 
TEST EVEN AFTER REDUCING H BY A FACTOR DGEA2190 

OF I.OEIO FROM ITS INITIAL VALUE. DGEA2200 

SEE REM7UUCS. DGEA2210 

lER = 133, IMPLIES THE INTEGRATION WAS DGEA2220 

HALTED AFTER FAILING TO ACHIEVE DGEA2230 

CORRECTOR CONVERGENCE EVEN AFTER DGEA2240 

REDUCING H BY A FACTOR OF I.OEIO FROM DGEA2250 
ITS INITIAL VALUE. SEE REMARKS. DGEA2260 

lER = 134, IMPLIES THAT AFTER SOME INITIAL DGEA2270 
SUCCESS, THE INTEGRATION WAS HALTED EITHERDGEA2280 
BY REPEATED ERROR TEST FAILURES OR BY DGEA2290 

A TEST ON TOL. SEE REMARKS. DGEA2300 

lER = 135, IMPLIES THAT ONE OF THE INPUT DGEA2310 
PARAMETERS N, X, H, XEND , TOL, METH, MITER, OR DGEA2320 
INDEX WAS SPECIFIED INCORRECTLY. DGEA2330 

lER = 13 6, IMPLIES THAT INDEX HAD A VALUE DGEA234 0 

OF -1 ON INPUT, BUT THE DESIRED CHANGES DGEA2350 
OF PARAMETERS WERE NOT IMPLEMENTED DGEA23 60 

BECAUSE XEND WAS NOT BEYOND X. DGEA2370 

INTERPOLATION TO X = XEND WAS PERFORMED. DGEA2380 
TO TRY AGAIN, SIMPLY CALL AGAIN WITH DGEA2390 

INDEX EQUAL TO -1 AND A NEW VALUE FOR DGEA2400 
XEND. DGEA2410 



Step - # of integration steps taken 

PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 

- SINGLE/H36,H48,H60 

REQD. IMSL ROUTINES - DGRCS , DGRIN, DGRPS , DGRST, LUDATF , LUELMF, LEQTIB , 



+ 

DGEA2420 

DGEA2430 

DGEA2440 

DGEA2450 

DGEA2460 
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NOTATION 



REMARKS 1 . 



2 . 



3 . 



UERTST , UGETIO DGEA2 47 0 

DGEA2480 

- INFORMATION ON SPECIAL NOTATION AND DGEA2490 

CONVENTIONS IS AVAILABLE IN THE MANUAL DGEA2500 

INTRODUCTION OR THROUGH IMSL ROUTINE UHELP DGEA2510 

DGEA2520 

THE EXTERNAL SUBROUTINE FCNJ IS USED ONLY WHEN DGEA2530 

INPUT PARAMETER MITER IS EQUAL TO 1 OR - 1 . OTHERWISE, DGEA2540 
A DUMMY FUNCTION CAN BE USED. THE DUMMY SUBROUTINE DGEA2550 

SHOULD BE OF THE FOLLOWING FORM DGEA25 60 

SUBROUTINE FCNJ (N,X,Y,PD) DGEA2570 

INTEGER N DGEA2580 

REAL Y(N) ,PD(N,N) ,X DGEA2590 

RETURN DGEA2600 

END DGEA2 610 

AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED DGEA2620 
(IER=0) AND A NORMAL CONTINUATION IS DESIRED, SIMPLY DGEA2630 
RESET XEND AND CALL DGEAR AGAIN. ALL OTHER DGEA264 0 

PARAMETERS WILL BE READY FOR THE NEXT CALL. A CHANGE DGEA2 650 

OF PARAMETERS WITH INDEX EQUAL TO -1 CAN BE MADE DGEA2660 

AFTER EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN. DGEA2 670 
THE COMMON BLOCKS /DBAND/ AND /GEAR/ NEED TO BE DGEA2680 

PRESERVED BETWEEN CALLS TO DGEAR. IF IT IS NECESSARY DGEA2690 

FOR THE COMMON BLOCKS TO EXIST IN THE CALLING PROGRAM DGEA2700 
THE FOLLOWING STATEMENTS SHOULD BE INCLUDED DGEA2710 

COMMON /DBAND/ NLC,NUC DGEA2720 

COMMON /GEAR/ DUMMY (48 ) , SDUMMY (4 ) , IDUMMY (3 8 ) DGEA2730 

WHERE DUMMY, SDUMMY, AND IDUMMY ARE VARIABLE NAMES NOT DGEA274 0 
USED ELSEWHERE IN THE CALLING PROGRAM. (FOR DOUBLE DGEA2750 
PRECISION DUMMY IS TYPE DOUBLE AND SDUMMY IS TYPE REAL) DGEA2760 
THE CHOICE OF VALUES FOR METH AND MITER MAY REQUIRE DGEA2770 
SOME EXPERIMENTATION, AND ALSO SOME CONSIDERATION OF DGEA27 80 

THE NATURE OF THE PROBLEM AND OF STORAGE REQUIREMENTS. DGEA2 790 
THE PRIME CONSIDERATION IS STIFFNESS. IF DGEA2800 

THE PROBLEM IS NOT STIFF, THE BEST CHOICE IS PROBABLY DGEA2810 
METH = 1 WITH MITER = 0 . IF THE PROBLEM IS STIFF TO A DGEA2820 
SIGNIFICANT DEGREE, THEN METH SHOULD BE 2 AND MITER DGEA2830 
SHOULD BE 1,2, -1,-2 OR 3 . IF THE USER HAS NO KNOWLEDGE DGEA2840 
OF THE INHERENT TIME CONSTANTS OF THE PROBLEM, WITH DGEA2850 
WHICH TO PREDICT ITS STIFFNESS, ONE WAY TO DETERMINE DGEA2860 
THIS IS TO TRY METH = 1 AND MITER = 0 FIRST, AND LOOK DGEA2 870 
AT THE BEHAVIOR OF THE SOLUTION COMPUTED AND THE STEP DGEA2880 
SIZES USED. IF THE TYPICAL VALUES OF H ARE MUCH DGEA2890 

SMALLER THAN THE SOLUTION BEHAVIOR WOULD SEEM TO DGEA2900 

REQUIRE (THAT IS, MORE THAN 100 STEPS ARE TAKEN OVER DGEA2910 
AN INTERVAL IN WHICH THE SOLUTIONS CHANGE BY LESS DGEA2920 

THAN ONE PERCENT) , THEN THE PROBLEM IS PROBABLY STIFF DGEA2930 
AND THE DEGREE OF STIFFNESS CAN BE ESTIMATED FROM THE DGEA294 0 
VALUES OF H USED AND THE SMOOTHNESS OF THE SOLUTION. DGEA2950 
IF THE DEGREE OF STIFFNESS IS ONLY SLIGHT, IT MAY BE DGEA2960 
THAT METH=1 IS MORE EFFICIENT THAN METH=2 . DGEA2970 

EXPERIMENTATION WOULD BE REQUIRED TO DETERMINE THIS. DGEA2980 
REGARDLESS OF METH, THE LEAST EFFECTIVE VALUE OF DGEA2990 

MITER IS 0, AND THE MOST EFFECTIVE IS 1,-1, 2, OR -2. DGEA3 000 
MITER = 3 IS GENERALLY SOMEWHERE IN BETWEEN. SINCE DGEA3 010 

THE STORAGE REQUIREMENTS GO UP IN THE SAME ORDER AS DGEA3 020 
EFFECTIVENESS, TRADE-OFF CONSIDERATIONS ARE DGEA3030 

NECESSARY. FOR REASONS OF ACCURACY AND SPEED, THE DGEA3 04 0 

CHOICE OF ABS (MITER) =1 IS GENERALLY PREFERRED TO DGEA3050 

ABS (MITER) =2 , UNLESS THE SYSTEM IS FAIRLY COMPLICATED DGEA3060 
(AND FCNJ IS THUS NOT FEASIBLE TO CODE) . THE DGEA3 070 

ACCURACY OF THE FCNJ CALCULATION CAN BE CHECKED BY DGEA3080 
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COMPARISON OF THE JACOBIAN WITH THAT GENERATED WITH 
ABS (MITER) =2 . IF THE JACOBIAN MATRIX IS SIGNIFICANTLY 
DIAGONALLY DOMINANT, THEN THE OPTION MITER = 3 IS 
LIKELY TO BE NEARLY AS EFFECTIVE AS ABS (MITER) =1 OR 2, 
AND WILL SAVE CONSIDERABLE STORAGE AND RUN TIME. 

IT IS POSSIBLE, AND POTENTIALLY QUITE DESIRABLE, TO 
USE DIFFERENT VALUES OF METH AND MITER IN DIFFERENT 
SUBINTERVALS OF THE PROBLEM. FOR EXAMPLE, IF THE 
PROBLEM IS NON- STIFF INITIALLY AND STIFF LATER, 

METH = 1 AND MITER = 0 MIGHT BE SET INITIALLY, AND 
METH = 2 AND MITER = 1 LATER. 

5. THE INITIAL VALUE OF THE STEP SIZE, H, SHOULD BE 
CHOSEN CONSIDERABLY SMALLER THAN THE AVERAGE VALUE 
EXPECTED FOR THE PROBLEM, AS THE FIRST-ORDER METHOD 
WITH WHICH DGEAR BEGINS IS NOT GENERALLY THE MOST 
EFFICIENT ONE. HOWEVER, FOR THE FIRST STEP, AS FOR 
EVERY STEP, DGEAR TESTS FOR THE POSSIBILITY THAT 
THE STEP SIZE WAS TOO LARGE TO PASS THE ERROR TEST 
(BASED ON TOL) , AND IF SO ADJUSTS THE STEP SIZE 
DOWN AUTOMATICALLY. THIS DOWNWARD ADJUSTMENT, IF 
ANY, IS NOTED BY lER HAVING THE VALUES GG OR 67, 

AND SUBSEQUENT RUNS ON THE SAME OR SIMILAR PROBLEM 
SHOULD BE STARTED WITH AN APPROPRIATELY SMALLER 
VALUE OF H. 

6 . SOME OF THE VALUES OF INTEREST LOCATED IN THE 

COMMON BLOCK /GEAR/ ARE 

A. PiUSED, THE STEP SIZE H LAST USED SUCCESSFULLY 

(DUMMY (8) ) 

B. NQUSED, THE ORDER LAST USED SUCCESSFULLY 

(IDUMMY(6) ) 

C. NSTEP, THE CUMULATIVE NUMBER OF STEPS TAKEN 

(IDUMMY(7) ) 

D. NFE, THE CUMULATIVE NUMBER OF FCN EVALUATIONS 

(IDUMMY(8) ) 

E. NJE, THE CUMULATIVE NUMBER OF JACOBIAN 
EVALUATIONS, AND HENCE ALSO OF MATRIX LU 
DECOMPOSITIONS (IDUMMYO)) 

7 . THE NORMAL USAGE OF DGEAR MAY BE SUMMARIZED AS FOLLOWS 

A. SET THE INITIAL VALUES IN Y. 

B. SET N, X, H, TOL, METH, AND MITER. 

C. SET XEND TO THE FIRST OUTPUT POINT, AND INDEX TO 1. 

D. CALL DGEAR 

E. EXIT IF lER IS GREATER THAN 128. 

F. OTHERWISE, DO DESIRED OUTPUT OF Y. 

G. EXIT IF THE PROBLEM IS FINISHED. 

H. OTHERWISE, RESET XEND TO THE NEXT OUTPUT POINT, AND 

RETURN TO STEP D. 

8 . THE ERROR WHICH IS CONTROLLED BY WAY OF THE PARAMETER 

TOL IS AN ESTIMATE OF THE LOCAL TRUNCATION ERROR, THAT 
IS, THE ERROR COMMITTED ON TAKING A SINGLE STEP WITH 
THE METHOD, STARTING WITH DATA REGARDED AS EXACT. THIS 
IS TO BE DISTINGUISHED FROM THE GLOBAL TRUNCATION 
ERROR, WHICH IS THE ERROR IN ANY GIVEN COMPUTED VALUE 
OF Y(X) AS A RESULT OF THE LOCAL TRUNCATION ERRORS 
FROM ALL STEPS TAKEN TO OBTAIN Y (X) . THE LATTER ERROR 
ACCUMULATES IN A NON -TRIVIAL WAY FROM THE LOCAL 
ERRORS, AND IS NEITHER ESTIMATED NOR CONTROLLED BY 
THE ROUTINE. SINCE IT IS USUALLY THE GLOBAL ERROR THAT 
A USER WANTS TO HAVE UNDER CONTROL, SOME 
EXPERIMENTATION MAY BE NECESSARY TO GET THE RIGHT 
VALUE OF TOL TO ACHIEVE THE USERS NEEDS . IF THE 
PROBLEM IS MATHEMATICALLY STABLE, AND THE METHOD USED 
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IS APPROPRIATELY STABLE, THEN THE GLOBAL ERROR AT A DGEA3710 
GIVEN X SHOULD VARY SMOOTHLY WITH TOL IN A MONOTONE DGEA3 720 
INCREASING MANNER. DGEA3730 

IF THE ROUTINE RETURNS WITH lER VALUES OF 132, 133, :'GEA3740 

OR 134, THE USER SHOULD CHECK TO SEE IF TOO MUCH DGEA3750 

ACCURACY IS BEING REQUIRED. THE USER MAY WISH TO DGEA3760 

SET TOL TO A LARGER VALUE AND CONTINUE. ANOTHER DGEA377 0 

POSSIBLE CAUSE OF THESE ERROR CONDITIONS IS AN DGEA3780 

ERROR IN THE CODING OF THE EXTERNAL FUNCTIONS FCN DGEA379 0 

OR FCNJ. IF NO ERRORS ARE FOUND, IT MAY BE NECESSARY DGEA3800 

TO MONITOR INTERMEDIATE QUANTITIES GENERATED BY THE DGEA3810 
ROUTINE. THESE QUANTITIES ARE STORED IN THE WORK VECTORDGEA3 820 
WK AND INDEXED BY SPECIFIC ELEMENTS IN THE COMMON BLOCKDGEA3 830 
/GEAR/. IF lER IS 132 OR 134, THE COMPONENTS CAUSING DGEA3840 
THE ERROR TEST FAILURE CAN BE IDENTIFIED FROM LARGE DGEA3850 
VALUES OF THE QUANTITY DGEA3860 

WK(IDUMMY(11) +1) /WK(I) , FOR 1=1 , . . . , N . DGEA3870 

ONE CAUSE OF THIS MAY BE A VERY SMALL BUT NONZERO DGEA3 880 

INITIAL VALUE OF ABS(Y(I)) . DGEA3890 

IF lER IS 133, SEVERAL POSSIBILITIES EXIST. DGEA3900 

IT MAY BE INSTRUCTIVE TO TRY DIFFERENT VALUES OF MITER . DGEA3 9 1 0 



COPYRIGHT 



WARRANTY 



ALTERNATIVELY, THE USER MIGHT MONITOR SUCCESSIVE 
CORRECTOR ITERATES CONTAINED IN WK (IDUMMY (12) +1) , FOR 
1=1, . . . ,N. ANOTHER POSSIBILITY MIGHT BE TO MONITOR 
THE JACOBIAN MATRIX, IF ONE IS USED, STORED, BY 
COLUMN, IN WK (IDUMMY (10) +1) , FOR 1=1 , . . . , N*N IF 
ABS (MITER) IS EQUAL TO 1 OR 2 , OR FOR 1 = 1, . . .,N IF 
MITER IS EQUAL TO 3. 

- 1984 BY IMSL, INC. ALL RIGHTS RESERVED. 
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DGEA3930 
DGEA3940 
DGEA3950 
DGEA3960 
DGEA3970 
DGEA3980 
DGEA3990 
DGEA4000 
DGEA4010 

IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN DGEA4020 
APPLIED TO THIS CODE. NO OTHER WARRANTY, DGEA4030 
EXPRESSED OR IMPLIED, IS APPLICABLE. DGEA4040 

DGEA4050 

• - DGEA4060 



DGEA4070 
DGEA4080 
+ 

DGEA4100 
+ 

DGEA4120 
DGEA4130 
DGEA4140 

JSTART , NSQ , NQUSED , NSTEP , NFE , NJE , I , NO , NHCUT , KGO , DGEA4 150 



SUBROUTINE DGEAR (N, FCN, FCNJ, X, H, Y, XEND , TOL, METH, MITER, INDEX, 
1 IWK,WK,IER,Step) 

SPECIFICATIONS FOR ARGUMENTS 
INTEGER N, METH, MITER, INDEX, I WK(1) , lER, step 

DOUBLE PRECISION X, H, Y (N) , XEND, TOL, WK ( 1) 

SPECIFICATIONS FOR LOCAL VARIABLES 
INTEGER NERROR , NSAVE 1 , NSAVE2 , NPW , NY , NC , MFC , KFLAG , 



JER,KER,NN,NEQUIL, IDUMMY (21) ,NLC,NUC 
REAL SDUMMY(4) 

DOUBLE PRECISION T, HH, HMIN, HMAX, EPSC , UROUND , EPS J, HUSED , TOUTP , 



EXTERNAL 
COMMON /DBAND/ 
COMMON /GEAR/ 



DGEA4160 

DGEA4170 

DGEA4180 

DGEA4190 

DGEA4200 

DGEA4210 

DGEA4220 

DGEA4230 



AYI,D,DN,SEPS, DUMMY (39) 

FCN, FCNJ 
NLC,NUC 

T , HH , HMIN , HMAX , EPSC , UROUND , EPS J , HUSED , DUMMY , 

1 TOUTP , SDUMMY , NC , MFC , KFLAG , JSTART , NSQ , NQUSED , 

2 NSTE P, NFE, NJE, NPW, NERROR, NSAVE 1, NSAVE 2, NEQUIL, DGEA4240 

3 NY , IDUMMY , NO , NHCUT DGEA4250 

DATA SEPS/Z3410000000000000/ DGEA4260 

FIRST EXECUTABLE STATEMENT DGEA4270 

IF (MITER. GE.O) NLC = -1 DGEA4280 

KER = 0 DGEA4290 

JER = 0 DGEA4300 

URCUND = SEPS DGEA4310 

COMPUTE WORK VECTOR INDICIES DGEA4320 
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on nnnonnoo on 



NERROR = N 
NSAVEl = NERROR+N 
NSAVE2 = NSAVEl +N 
NY = NSAVE2+N 

IF (METH.EQ.l) NEQUIL = m + 13*N 
IF (METH.EQ.2) NEQUIL = NY+6*N 
NPW = NEQUIL + N 

IF (MITER. EQ. 0 .OR. MITER. EQ. 3) NPW = NEQUIL 
MFC = 10*METH+IABS (MITER) 

CHECK FOR INCORRECT INPUT PARAMETERS 



IF 

IF 

IF 

IF 

IF 

IF 

IF 

IF 

IF 

IF 



(MITER.lt. -2 .OR. MITER. GT. 3) GOTO 85 

(METH.NE.l .AND.METH.NE.2) GO TO 85 

(TOL.LE.O.DO) GO TO 85 

(N.LE.O) GO TO 85 

( (X-XEND) *H.GE.0.D0) GO TO 85 

(INDEX. EQ.O) GO TO 10 

(INDEX. EQ. 2) GO TO 15 

(INDEX. EQ. -1) GO TO 20 

(INDEX.EQ.3) GO TO 25 

( INDEX. NE.l) GO TO 85 



IF INITIAL VALUES OF YMAX OTHER THAN 
THOSE SET BELOW ARE DESIRED, THEY 
SHOULD BE SET HERE. ALL YMAX (I) 
MUST BE POSITIVE. IF VALUES FOR 
HMIN OR HMAX, THE BOUNDS ON 
DABS (HH) , OTHER THAN THOSE BELOW 
ARE DESIRED, THEY SHOULD BE SET 
BELOW . 



DO 



5 1=1, N 
WK(I) = DABS (Y(I) ) 

IF (WK(I) .EQ.O.DO) WK(I) = l.DO 
WK(NY+I) = Y(I) 

CONTINUE 
NC = N 
T = X 
HH = H 

IF ( (T+HH) .EQ.T) KER = 33 
HMIN = DABS(H) 

HMAX = DABS (X-XEND) *10 .DO 

EPSC = TOL 

JSTART = 0 

NO = N 

NSQ = N0*N0 

EPSJ = DSQRT (UROUND) 

NHCUT = 0 
DUMMY(2) = l.ODO 
DUMMY (14) = l.ODO 

GO TO 30 



10 HMAX 

GO TO 45 



TOUTP IS THE PREVIOUS VALUE OF XEND 
FOR USE IN HMAX. 

DABS (XEND -TOUTP) *10 .DO 



15 



20 



HMAX = DABS (XEND -TOUTP) *10 .DO 
IF ( (T-XEND) *HH.GE.0.D0) GO TO 95 
GO TO 50 



IF ( (T-XEND) *HH.GE.0. DO) 
JSTART = -1 
NC = N 
EPSC = TOL 



GO TO 90 
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25 IF ( (T+HH) .EQ.T) KER = 33 

write {*,*),' error code = ' , ker 

30 NN = NO 

step = step + 1 

write {*,*)' step = '.step 

CALL DGRST {FCN, FCNJ, WK (NY+1) , WK, WK (NERROR+1 ) , WK (NSAVEl+1) , 
1 WK{NSAVE2+1) ,WK{NPW+1) , WK (NEQUIL+1) , IWK, NN, step) 



KGO = 1-KFLAG 
GO TO {35,55,70,80) 

35 CONTINUE 



KGO 



KFLAG =0, -1, -2, -3 



C NORMAL RETURN FROM INTEGRATOR. THE 

C WEIGHTS YMAX(I) ARE UPDATED. IF 

C DIFFERENT VALUES ARE DESIRED, THEY 

C SHOULD BE SET HERE. A TEST IS MADE 

C FOR TOL BEING TOO SMALL FOR THE 

C MACHINE PRECISION. ANY OTHER TESTS 

C OR CALCULATIONS THAT ARE REQUIRED 

C AFTER EVERY STEP SHOULD BE 

C INSERTED HERE. IF INDEX =3, Y IS 

C SET TO THE CURRENT SOLUTION ON 

C RETURN. IF INDEX =2, HH IS 

C CONTROLLED TO HIT XEND (WITHIN 

C ROUNDOFF ERROR) , AND THEN THE 

C CURRENT SOLUTION IS PUT IN Y ON 

C RETURN. FOR ANY OTHER VALUE OF 

C INDEX, CONTROL RETURNS TO THE 

C INTEGRATOR UNLESS XEND HAS BEEN 

C REACHED. THEN INTERPOLATED VALUES 

C OF THE SOLUTION ARE COMPUTED AND 

C STORED IN Y ON RETURN. 

C IF INTERPOLATION IS NOT 

C DESIRED, THE CALL TO DGRIN SHOULD 

C BE REMOVED AND CONTROL TRANSFERRED 

C TO STATEMENT 95 INSTEAD OF 105. 

D = O.DO 
DO 40 1=1, N 

AYI = DABS (WK{NY+I) ) 

WK(I) = DMAXKWK(I) ,AYI) 

40 D = D+ (AYI/WK(I) ) **2 
D = D* (UROUND/TOL) **2 
DN = N 

IF (D.GT.DN) GO TO 75 
IF (INDEX. EQ. 3) GO TO 95 
IF (INDEX. EQ. 2) GO TO 50 
45 IF ( (T-XEND) *HH.LT. 0 .DO) GO TO 25 
NN = NO 

CALL DGRIN (XEND, WK (NY+1) ,NN,Y) 

X = XEND 
GO TO 105 

50 IF ({ (T+HH) -XEND) *HH.LE. O.DO) GO TO 25 

IF (DABS (T-XEND) .LE.UROUND*DMAXl (10. D0*DABS (T) ,HMAX) ) GO TO 95 
IF { (T-XEND) *HH.GE. O.DO) GO TO 95 
HH = (XEND-T) * (1 .DO-4 .D0*UROUND) 

JSTART = -1 
GO TO 25 

C ON AN ERROR RETURN FROM INTEGRATOR, 

C AN IMMEDIATE RETURN OCCURS IF 

C KFLAG = -2, AND RECOVERY ATTEMPTS 
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65 
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JER = 66 

IF (NHCUT.EQ. 10) GO TO 65 

NHCUT = NHCUT+1 

HMIN = HMIN* . IDO 

HH = HH*.1D0 

JSTART = -1 

GO TO 25 

IF (JER.EQ.66) JER = 132 
IF (JER.EQ.67) JER = 133 
GO TO 95 

JER = 134 
GO TO 95 



ARE MADE OTHERWISE. TO RECOVER HH 
and HMIN ARE REDUCED BY A FACTOR 
. 1 UP TO 10 TIMES BEFORE GIVING 



75 



80 



85 



90 



JER = 134 
KFLAG = -2 
GO TO 95 

JER = 67 
GO TO 60 

JER = 135 
GO TO 110 

JER = 136 
NN = NO 

CALL DGRIN (XEND, WK (NY+1) NN Y) 
X = XEND ' ' 

GO TO 110 



95 



X = T 

DO 100 1=1, N 
Y(I) = WK(NY+I) 

IF (JER. LT. 128) 

TOUTP = X 
IF (KFLAG. EQ.O) 

IF (KFLAG. NE.O) .. _ 
lER = MAXO (KER, JER) 

CONTINUE 

p (KER. NE.O. AND. JER. LT. 128) CALL UE RTS T (KER fiHnruno \ 
9005 UERTST ( JER, 6HDGEAR ) 

END 



100 

105 



110 

9000 



INDEX = KFLAG 

H = HUSED 
H = HH 
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IMSL ROUTINE NAME - DGRST DGRSOOlO 

DGRS0020 

•modified to print sheath and diagnostic output to files ’’sheatha.dat” + 
and ’’diag.dat" 

COMPUTER 



LATEST REVISION 



PURPOSE 

PRECISION/HARDWARE 



REQD. IMSL ROUTINES - 



NOTATION 



- IBM/DOUBLE 

- JUNE 1, 1982 

- NUCLEUS CALLED ONLY BY IMSL SUBROUTINE DGEAR 



SINGLE AND DOUBLE /H3 2 
SINGLE/H36,H48,H60 

DGRCS , DGRPS , LUDATF , LUELMF , LEQTIB , UERTST , 
UGETIO 

INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 



COPYRIGHT 



WARRANTY 



- 1982 BY IMSL, INC. ALL RIGHTS RESERVED. 



SUBROUTINE DGRST 



INTEGER 

DOUBLE PRECISION 



INTEGER 



1 

2 

3 

4 



REAL 

DOUBLE PRECISION 



1 

2 

3 



EXTERNAL 
COMMON /DBAND/ 
COMMON /GEAR/ 



1 

2 

3 

4 

5 

6 

7 

8 



open (unit=8 , file 
open (unit=9 , file 
KFLAG = 0 
TOr^D = T 



DGRS0050 
DGRS0060 
DGRS0070 
DGRS0080 
DGRS0090 
DGRSOlOO 
DGRSOllO 
DGRS0120 
DGRS0130 
DGRS0140 
DGRS0150 
DGRS0160 
DGRS0170 
DGRS0180 
DGRS0190 
DGRS0200 
DGRS0210 
DGRS0220 

IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN DGRS0230 
APPLIED TO THIS CODE. NO OTHER WARRANTY, DGRS0240 

EXPRESSED OR IMPLIED, IS APPLICABLE. DGRS0250 

DGRS0260 

DGRS0270 

DGRS0280 
DGRS0290 
+ 

DGRS0310 
DGRS0320 
DGRS0330 
+ 

DGRS0350 
DGRS0360 
DGRS0370 
DGRS0380 
DGRS0390 
DGRS0400 
DGRS0410 
DGRS0420 
DGRS0430 
,DGRS0440 
+ 

DGRS0460 
DGRS0470 
DGRS0480 
DGRS0490 
DGRS0500 
DGRS0510 
DGRS0520 
DGRS0530 
DGRS0540 
DGRS0550 
DGRS0560 
DGRS0570 
+ 

+ 

DGRS0580 
DGRS0590 
DGRS0600 



(FCN, FCNJ, Y, YMAX, ERROR, SAVEl , SAVE2 , PW, EQUIL, 
IPIV,NO,Step) 

SPECIFICATIONS FOR ARGUMENTS 
IPIV(l) ,N0 

Y(N0,1) ,YMAX(1) , ERROR (1) , SAVEl (1) ,SAVE2(1) , 

PW ( 1 ) , EQUIL ( 1 ) , epr ime , epr ime ( 2 ) 

SPECIFICATIONS FOR LOCAL VARIABLES 
N , MF , KFLAG , JSTART , NQUSED , NSTEP , NFE , NJE , NSQ , 

I , METH , MITER , NQ , L , IDOUB , MFOLD , NOLD , IRET , MEO , 
MIO , I WEVAL , MAXDER , LMAX , IREDO , J , NSTEP J , J1 , J2 , 

M, lER, NEWQ, NPW, NERROR, NSAVEl , NSAVE2 , NEQUIL, NY, 
MITERl , IDUMMY ( 2 ) , NLC , NUC , NWK , JER 

TQ(4) 

T,H,HMIN,HMAX,EPS,UR0UND,HUSED,EL(13) , OLDLO , 
TOLD , RMAX , RC , CRATE , EPSOLD , HOLD , FN , EDN , E , EUP , 
BND , RH, R1 , CON, R, HLO , RO , D , PHLO , PR3 , D1 , ENQ3 , ENQ2 
PR2 , PRl , ENQl , EPS J , DUMMY , t cum 
FCN, FCNJ 
NLC , NUC 

T , H , HMIN , HMAX , EPS , UROUND , EPS J , HUSED , 

EL , OLDLO , TOLD , RMAX , RC , CRATE , EPSOLD , HOLD , FN , 
EDN, E , EUP , BND , RH, R1 , R, HLO , RO , D , PHLO , PR3 , D1 , 
ENQ3 , ENQ2 , PR2 , PRl , ENQl , DUMMY , TQ , 

N , MF , KFLAG , JSTART , NSQ , NQUSED , NSTE P , NFE , NJE , 
NPW , NERROR , NSAVEl , NSAVE2 , NEQUIL , NY , 

I , METH , MITER , NQ , L , IDOUB , MFOLD , NOLD , IRET , MEO , 
MIO , IWEVAL , MAXDER , LMAX , IREDO , J , NSTEP J , J1 , J2 , 
M,NEWQ, IDUMMY 

FIRST EXECUTABLE STATEMENT 
=' sheatha.dat' , status= ' unknown' ) 

=' diag.dat' , status= ' unknown' ) 



THIS ROUTINE PERFORMS ONE STEP OF 
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nooononoooo non 



IF (JSTART.GT.O) GO TO 50 
IF (JSTART.NE. 0) GO TO 10 



THE INTEGRATION OF AN INITIAL 
VALUE PROBLEM FOR A SYSTEM OF 
ORDINARY DIFFERENTIAL EQUATIONS. 



C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



ON THE FIRST CALL, THE ORDER IS SET 
TO 1 AND THE INITIAL YDOT IS 
CALCULATED. RMAX IS THE MAXIMUM 
RATIO BY WHICH H CAN BE INCREASED 
IN A SINGLE STEP. IT IS INITIALLY 
1.E4 TO COMPENSATE FOR THE SMALL 
INITIAL H, BUT THEN IS NORMALLY 
EQUAL TO 10. IF A FAILURE OCCURS 
(IN CORRECTOR CONVERGENCE OR ERROR 
TEST) , RMAX IS SET AT 2 FOR THE 
NEXT INCREASE. 

CALL FCN (N, T, Y, SAVEl , eprime, eprime2) 

DO 5 1=1, N 

5 Y (I, 2) = H*SAVE1 (I) 

METH = MF/10 

MITER = MF-10*METH 

NQ = 1 

L = 2 

IDOUB = 3 

RMAX = 1.D4 

RC = O.DO 

CRATE = l.DO 

HOLD = H 

MFOLD = MF 

NSTEP = 0 

NSTEPJ = 0 

NFE = 1 

NJE = 0 

IRET = 3 

GO TO 15 

IF THE CALLER HAS CHANGED METH, 

DGRCS IS CALLED TO SET THE 
COEFFICIENTS OF THE METHOD. IF THE 
CALLER HAS CHANGED N, EPS, OR 
METH, THE CONSTANTS E, EDN, EUP, 
AND BND MUST BE RESET. E IS A 
COMPARISON FOR ERRORS OF THE 
CURRENT ORDER NQ. EUP IS TO TEST 
FOR INCREASING THE ORDER, EDN FOR 
DECREASING THE ORDER. BND IS USED 
TO TEST FOR CONVERGENCE OF THE 
CORRECTOR ITERATES. IF THE CALLER 
HAS CHANGED H, Y MUST BE RESCALED. 
IF H OR METH HAS BEEN CHANGED, 
IDOUB IS RESET TO L + 1 TO PREVENT 
FURTHER CHANGES IN H FOR THAT MANY 
STEPS . 

10 IF (MF.EQ. MFOLD) GO TO 25 
MEO = METH 
MIO = MITER 
METH = MF/10 
MITER = MF-10*METH 
MFOLD = MF 

IF (MITER. NE. MIO) IWEVAL = MITER 
IF (METH. EQ. MEO) GO TO 25 
IDOUB = L+1 
IRET = 1 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



25 



30 



35 

40 



45 



50 



55 



C 

C 

C 

C 

C 

C 

C 

C 

C 



CALL DGRCS (METH , NQ , EL, TQ , MAXDER) 

LMAX = MAXDER+1 
RC = RC*EL(1) /OLDLO 
OLDLO = EL(1) 

FN = N 

EDN = FN* (TQ (1) *EP£) *^^2 
E = FN* (TQ (2) *EPS) **2 
EUP = FN* (TQ{3) *EPS) **2 
END = FN* (TQ(4) *EPS) **2 
EPSOLD = EPS 
NOLD = N 

GO TO (30,35,50) , IRET 

IF ( (EPS. EQ. EPSOLD) .AND. (N.EQ.NOLD)) GO TO 30 
IF (N.EQ.NOLD) IWEVAL = MITER 
IRET = 1 
GO TO 20 

IF (H.EQ.HOLD) GO TO 50 
RH = H/HOLD 
H = HOLD 
IREDO = 3 
GO TO 40 

RH = DMAXl (RH,HMIN/DABS (H) ) 

RH = DMINl (RH,HMAX/DABS (H) ,RMAX) 

R1 = l.DO 
DO 45 J=2,L 
R1 = R1*RH 
DO 45 1=1, N 
Y(I,J) = Y(I,J)*R1 
H = H*RH 
RC = RC*RH 
IDOUB = L+1 

IF (IREDO. EQ.O) GO TO 285 

THIS SECTION COMPUTES THE PREDICTED 
VALUES BY EFFECTIVELY MULTIPLYING 
THE Y ARRAY BY THE PASCAL TRIANGLE 
MATRIX. RC IS THE RATIO OF NEW TO 
OLD VALUES OF THE COEFFICIENT 
H*EL(1) . WHEN RC DIFFERS FROM 1 BY 
MORE THAN 30 PERCENT, OR THE 
CALLER HAS CHANGED MITER, IWEVAL 
IS SET TO MITER TO FORCE THE 
PARTIALS TO BE UPDATED, IF 
PARTIALS ARE USED. IN ANY CASE, 

THE PARTIALS ARE UPDATED AT LEAST 
EVERY 20 -TH STEP. 

IF (DABS (RC-1 .DO) .GT.0.3D0) IWEVAL = MITER 
IF (NSTEP.GE.NSTEPJ+20) IWEVAL = MITER 
T = T+H 
DO 55 J1=1,NQ 
DO 55 J2=J1,NQ 
J = (NQ+Jl) -J2 
DO 55 1=1, N 

Y(I,J) = Y(I, J)+Y(I, J+1) 

UP TO 3 CORRECTOR ITERATIONS ARE 
TAKEN. A CONVERGENCE TEST IS MADE 
ON THE R.M.S. NORM OF EACH 
CORRECTION, USING BND, WHICH IS 
DEPENDENT ON EPS. THE SUM OF THE 
CORRECTIONS IS ACCUMULATED IN THE 
VECTOR ERROR (I) . THE Y ARRAY IS 
NOT ALTERED IN THE CORRECTOR LOOP. 
THE UPDATED Y VECTOR IS STORED 
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C TEMPORARILY IN SAVEl . DGRS1850 

60 DO 65 1=1, N DGRS1860 

65 ERROR (I) = O.DO DGRS1870 

M = 0 DGRS1880 

CALL FCN (N,T, Y, SAVE2 , eprime, eprime2) + 

NFE = NFE+1 DGRS1900 

IF (IWEVAL.LE.O) GO TO 95 DGRS1910 

C IF INDICATED, THE MATRIX P = I - DGRS1920 

C H*EL(1)*J IS REEVALUATED BEFORE DGRS1930 

C STARTING THE CORRECTOR ITERATION. DGRS1940 

C IWEVAL IS SET TO 0 AS AN INDICATOR DGRS1950 

C THAT THIS HAS BEEN DONE. IF MITER DGRS1960 

C = 1 OR 2, P IS COMPUTED AND DGRS1970 

C PROCESSED IN PSET. IF MITER = 3, DGRS1980 

C THE MATRIX USED IS P = I - DGRS1990 

C H*EL(1)*D, WHERE D IS A DIAGONAL DGRS2000 

C MATRIX. DGRS2010 

IWEVAL = 0 DGRS2020 

RC = l.DO DGRS2030 

NJE = NJE+1 DGRS2040 

NSTEPJ = NSTEP DGRS2050 

GOTO (75,70,80), MITER DGRS2060 

70 NFE = NFE+N DGRS2070 

75 CON = -H*EL(1) DGRS2080 

MITERl = MITER DGRS2090 

CALL DGRPS (FCN, FCNJ, Y, NO , CON, MITERl , YMAX, SAVEl , SAVE2 , PW, EQUIL, DGRS2100 

1 IPIV,IER) DGRS2110 

IF (lER.NE.O) GO TO 155 DGRS2120 

GO TO 125 DGRS2130 

80 R = EL (1)*. IDO DGRS2140 

DO 85 1=1, N DGRS2150 

85 PW(I) = Y(I, 1) +R* (H*SAVE2 (I) -Y(I,2) ) DGRS2160 

CALL FCN (N, T, PW, SAVEl, eprime, eprime2) + 

NFE = NFE+1 DGRS2180 

HLO = H*EL(1) DGRS2190 

DO 90 1=1, N DGRS2200 

RO = H*SAVE2 (I) -Y(I,2) DGRS2210 

PW(I) = l.DO DGRS2220 

D= .1D0*R0-H* (SAVEl (I) -SAVE2 (I) ) DGRS2230 

SAVEl (I) = O.DO DGRS2240 

IF (DABS (RO) .LT.UROUND*YMAX(I) ) GO TO 90 DGRS2250 

IF (DABS (D) .EQ. O.DO) GOTO 155 DGRS2260 

PW(I) = .1D0*R0/D DGRS2270 

SAVEl (I) = PW(I)*R0 DGRS2280 

90 CONTINUE DGRS2290 

GO TO 135 DGRS2300 

95 IF (MITER. NE.O) GOTO (125,125,105), MITER DGRS2310 

DGRS2320 

IN THE CASE OF FUNCTIONAL ITERATION, DGRS2330 

UPDATE Y DIRECTLY FROM THE RESULT DGRS2340 

OF THE LAST FCN CALL. DGRS2350 

D = O.DO DGRS2360 

DO 100 1=1, N DGRS2370 

R = H*SAVE2 (I) -Y(I,2) DGRS2380 

D = D+ ( (R-ERROR(I) ) /YMAX(I) ) **2 DGRS2390 

SAVEl(I) = Y (I, 1) +EL (1) *R DGRS2400 

100 ERROR (I) = R DGRS2410 

GO TO 145 DGRS2420 

IN THE CASE OF THE CHORD METHOD, DGRS2430 

COMPUTE THE CORRECTOR ERROR, F SUB DGRS2440 
(M) , AND SOLVE THE LINEAR SYSTEM DGRS2450 

WITH THAT AS RIGHT-HAND SIDE AND P DGRS2460 
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105 



110 

115 

120 

125 

130 



131 

135 



140 



145 



AS COEFFICIENT MATRIX, USING THE 
LU DECOMPOSITION IF MITER = 1 OR 
2. IF MITER = 3, THE COEFFICIENT 
H*EL{1) IN P IS UPDATED. 

PHLO = HLO 
HLO = H*EL(1) 

IF {HLO .EQ. PHLO) GO TO 115 
R = HLO /PHLO 
DO 110 1 = 1, N 

D = 1 .DO-R* (1 .DO-1 .D0/PW(I) ) 

IF (DABS(D) .EQ.O.DO) GO TO 165 
PW(I) = l.DO/D 
DO 120 1=1, N 

SAVEl(I) = PW(I) * (H*SAVE2 (I) - (Y(I, 2) +ERROR(I) ) ) 

GO TO 135 
DO 130 1 = 1, N 

SAVEl(I) = H*SAVE2 (I) - (Y(I,2) +ERROR(I) ) 

IF (NLC .EQ. -1) GO TO 131 
NWK = (NLC+NUC+1) *N0+1 

CALL LEQT1B(PW,N,NLC,NUC,N0,SAVE1,1,N0,2,PW(NWK) , JER) 

GO TO 135 

CALL LUELMF (PW, SAVEl , IPIV, N, NO , SAVEl ) 

D = O.DO 
DO 140 1 = 1, N 

ERROR (I) = ERROR(I) +SAVEKD 

D = D+ (SAVEl(I) /YMAX(I) ) **2 
SAVEl(I) = Y(I,1) +EL(1) *ERROR(I) 

TEST FOR CONVERGENCE. IF M.GT.O, THE 
SQUARE OF THE CONVERGENCE RATE 
CONSTANT IS ESTIMATED AS CRATE, 

AND THIS IS USED IN THE TEST. 

IF (M.NE.O) CRATE = DMAXl ( . SDO’^CRATE , D/Dl ) 

IF ( (D*DMIN1(1 .DO, 2 .D0*CRATE) ) .LE.BND) GO TO 170 
D1 = D 
M = M+1 

IF (M.EQ.3) GO TO 150 

CALL FCN (N,T, SAVEl, SAVE2,eprime,eprime2) 

GO TO 95 

THE CORRECTOR ITERATION FAILED TO 
CONVERGE IN 3 TRIES. IF PARTIALS 
ARE INVOLVED BUT ARE NOT UP TO 
DATE, THEY ARE REEVALUATED FOR THE 
NEXT TRY. OTHERWISE THE Y ARRAY IS 
RETRACTED TO ITS VALUES BEFORE 
PREDICTION, AND H IS REDUCED, IF 
POSSIBLE. IF NOT, A NO- CONVERGENCE 
EXIT IS TAKEN. 



150 

155 



160 



165 



NFE = NFE+2 

IF (IWEVAL.EQ. -1) GO TO 165 
T = TOLD 
RMAX = 2. DO 
DO 160 J1=1,NQ 
DO 160 J2=J1,NQ 
J = (NQ+JD-J2 
DO 160 1=1, N 
Y(I, J) = Y(I, J) -Y(I, J+1) 

IF (DABS (H) .LE .HMIN*1 .OOOOIDO) 

RH = .25D0 

I REDO = 1 

GO TO 35 

IWFVAL = MITER 

GO ‘TO 60 



GO TO 280 
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THE CORRECTOR HAS CONVERGED. IWEVAL 
IS SET TO -1 IF PARTIAL 
DERIVATIVES WERE USED, TO SIGNAL 
THAT THEY MAY NEED UPDATING ON 
SUBSEQUENT STEPS. THE ERROR TEST 
IS MADE AND CONTROL PASSES TO 
STATEMENT 190 IF IT FAILS. 



170 IF (MITER. NE.O) IWEVAL = -1 
NFE = NFE+M 
D = O.DO 
DO 175 1 = 1, N 

175 D = D+ (ERROR(I) /YMAX(I) ) 

IF (D.GT.E) GO TO 190 



C 

C 

C 

C 

C 

C 

C 



AFTER A SUCCESSFUL STEP, UPDATE THE 
Y ARRAY. CONSIDER CHANGING H IF 
IDOUB = 1. OTHERWISE DECREASE 
IDOUB BY 1. IF IDOUB' IS THEN 1 AND 
NQ .LT. MAXDER, THEN ERROR IS 
SAVED FOR USE IN A POSSIBLE ORDER 
INCREASE ON THE NEXT STEP. IF A 
CHANGE IN H IS CONSIDERED, AN 
INCREASE OR DECREASE IN ORDER BY 
ONE IS CONSIDERED ALSO. A CHANGE 
IN H IS MADE ONLY IF IT IS BY A 
FACTOR OF AT LEAST 1.1. IF NOT, 
IDOUB IS SET TO 10 TO PREVENT 
TESTING FOR THAT MANY STEPS. 



KFLAG = 0 
IREDO = 0 
NSTEP = NSTEP+1 
HUSED = H 
NQUSED = NQ 
DO 180 J=1,L 
DO 180 1=1, N 

180 Y(I,J) = Y(I, J) +EL(J) *ERROR(I) 
IF ( IDOUB. EQ.l) GO TO 200 
IDOUB = IDOUB -1 
IF (IDOUB. GT.l) GO TO 290 
IF (L.EQ.LMAX) GO TO 290 
DO 185 1=1, N 

185 Y(I,LMAX) = ERROR(I) 

GO TO 290 



THE ERROR TEST FAILED. KFLAG KEEPS 
TRACK OF MULTIPLE FAILURES. 

RESTORE T AND THE Y ARRAY TO THEIR 
PREVIOUS VALUES, AND PREPARE TO 
TRY THE STEP AGAIN. COMPUTE THE 
OPTIMUM STEP SIZE FOR THIS OR ONE 
LOWER ORDER. 



190 



DO 195 
J = 
DO 195 
195 Y(I,J) 
RMAX = 
IF 
IF 



KFLAG = KFLAG -1 
T = TOLD 
DO 195 J1=1,NQ 
J2=J1,NQ 
(NQ+Jl) -J2 
1 = 1, N 

= Y(I, J) -Y(I, J+1) 

2. DO 

(DABS(H) .LE.HMIN^l.OOOOlDO) 
(KFLAG.LE. -3) GO TO 260 
IREDO = 2 
PR3 = l.D+20 
GO TO 210 



GO TO 270 
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REGARDLESS OF THE SUCCESS OR FAILURE 
OF THE STEP, FACTORS PRl , PR2 , AND 
PR3 ARE COMPUTED, BY WHICH H COULD 
BE DIVIDED AT ORDER NQ - 1 , ORDER 
NQ, OR ORDER NQ + 1 , RESPECTIVELY. 
IN THE CASE OF FAILURE, PR3 = 

1.E20 TO AVOID AN ORDER INCREASE. 
THE SMALLEST OF THESE IS 
DETERMINED AND THE NEW ORDER 
CHOSEN ACCORDINGLY. IF THE ORDER 
IS TO BE INCREASED, WE COMPUTE ONE 
ADDITIONAL SCALED DERIVATIVE. 



200 PR3 = 1 .D+20 

IF (L.EQ.LMAX) GO TO 210 
D1 = O.DO 
DO 205 1=1, N 

205 D1 = D1+ ( (ERROR(I) -Y(I,LMAX) ) /YMAX(I) ) **2 
ENQ3 = .5D0/(L+1) 

PR3 = ( (Dl/EUP) **ENQ3) *1 .4D0+1 .4D-6 
210 ENQ2 = .5D0/L 

PR2 = ( (D/E) **ENQ2) *1.2D0+1.2D-6 
PRl = l.D+20 
IF (NQ.EQ.l) GO TO 220 
D = O.DO 
DO 215 1=1, N 

215 D = D+ (Y(I,L) /YMAX(I) ) **2 
ENQl = .5D0/NQ 

PRl = ( (D/EDN) **ENQ1) *1 .3D0+1 .3D-6 
220 IF (PR2.LE.PR3) GO TO 225 
IF (PR3 .lt. PR l) GO TO 235 
GO TO 230 

225 IF (PR2.GT.PR1) GO TO 230 
NEWQ = NQ 
RH = 1.D0/PR2 
GO TO 250 
230 NEWQ = NQ-1 
RH = l.DO/PRl 

IF (KFLAG.NE.O.AND.RH.GT.l.DO) RH = l.DO 
GO TO 250 
235 NEWQ = L 

RH = 1.D0/PR3 

IF (RH.LT.l.lDO) GO TO 245 
DO 240 1 = 1, N 

240 Y(I,NEWQ+1) = ERROR ( I ) *EL (L) /L 
GO TO 255 
245 IDOUB = 10 
GO TO 290 

250 IF ( (KFLAG.EQ.O) .AND. (RH.LT.l .IDO) ) GO TO 245 



IF (NEWQ.EQ.NQ) GO TO 35 
255 NQ = NEWQ 
L = NQ+1 
IRET = 2 
GO TO 15 



IF THERE IS A CHANGE OF ORDER, RESET 
NQ, L, AND THE COEFFICIENTS. IN 
ANY CASE H IS RESET ACCORDING TO 
RH AND THE Y ARRAY IS RESCALED. 
THEN EXIT FROM 285 IF THE STEP WAS 
OK, OR REDO THE STEP OTHERWISE. 



CONTROL REACHES THIS SECTION IF 3 OR 
MORE FAILURES HAVE OCCURED . IT IS 
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ASSUMED THAT THE DERIVATIVES THAT 
HAVE ACCUMULATED IN THE Y ARRAY 
HAVE ERRORS OF THE WRONG ORDER. 
HENCE THE FIRST DERIVATIVE IS 
RECOMPUTED, AND THE ORDER IS SET 
TO 1. THEN H IS REDUCED BY A 
FACTOR OF 10, AND THE STEP IS 
RETRIED. AFTER A TOTAL OF 7 
FAILURES, AN EXIT IS TAKEN WITH 
KFLAG = -2. 

260 IF (KFLAG. EQ. -7) GO TO 275 
RH = .IDO 

RH = DMAXl (HMIN/DABS (H) ,RH) 

H = H*RH 

CALL FCN (N,T, Y,SAVEl,eprime,eprime2) 

NFE = NFE+1 
DO 265 1=1, N 
265 Y(I,2) = H*SA\rEl (I) 

IWEVAL = MITER 
IDOUB =10 

IF (NQ.EQ.l) GO TO 50 
NQ = 1 
L = 2 
IRET = 3 
GO TO 15 

ALL RETURNS ARE MADE THROUGH THIS 
SECTION. H IS SAVED IN HOLD TO 
ALLOW THE CALLER TO CHANGE H ON 
THE NEXT STEP. 

270 KFLAG = -1 
GO TO 290 
275 KFLAG = -2 
GO TO 290 
280 KFLAG = -3 
GO TO 290 
285 RMAX = 10. DO 
290 HOLD = H 

JSTART = NQ 

C- -Diagnostic Check of first and second derivatives of E 
if (team. eq. told) go to 310 

write (8,300) tcum, step,y (1, 1) ,y(2,l) ,y(3,l) ,y(4,l) ,y(5,l) 

300 format ( lx, ell. 4, lx, 15, 5( lx, ell. 4) ) 
write (9 , 305) step, eprime , eprime2 
305 format(lx,I5,2 (lx,e20.13) ) 

RETURN 

END 



DGRS4330 

DGRS4340 

DGRS4350 

DGRS4360 

DGRS4370 

DGRS4380 

DGRS4390 

DGRS4400 

DGRS4410 

DGRS4420 

DGRS4430 

DGRS4440 

DGRS4450 

DGRS4460 

+ 

DGRS4480 

DGRS4490 

DGRS4500 

DGRS4510 

DGRS4520 

DGRS4530 

DGRS4540 

DGRS4550 

DGRS4560 

DGRS4570 

DGRS4580 

DGRS4590 

DGRS4600 

DGRS4610 

DGRS4620 

DGRS4630 

DGRS4640 

DGRS4650 

DGRS4660 

DGRS4670 

DGRS4680 

DGRS4690 

DGRS4700 

+ 

+ 

+ 

+ 

+ 

+ 

DGRS4710 

DGRS4720 
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c 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 



IMSL ROUTINE NAME 



- DGRCS 



COMPUTER 
LATEST REVISION 
PURPOSE 

PRECISION/HARDWARE 

REQD . IMSL ROUTINES 
NOTATION 



- IBM/DOUBLE 

- JANUARY 1, 1978 

- NUCLEUS CALLED ONLY BY IMSL SUBROUTINE DGEAR 

- SINGLE AND DOUBLE /H3 2 

- SINGLE/H36,H48,H60 



COPYRIGHT 

WARRANTY 



- NONE REQUIRED 

- INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 

- 1978 BY IMSL, INC. ALL RIGHTS RESERVED. 

- IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 

APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 



SUBROUTINE DGRCS 



INTEGER 

REAL 

DOUBLE PRECISION 

INTEGER 

REAL 

DATA 



(METH , NQ , EL , TQ , MAXDER) 

SPECIFICATIONS FOR ARGUMENTS 
METH, NQ, MAXDER 
TQ(1) 

EL(1) 

SPECIFICATIONS FOR LOCAL VARIABLES 
K 

PERTST(12,2,3) 

PERTST/1 . , 1 . , 2 . , 1 . , .3158 , . 7407E-1, 

1 .1391E-1, .2182E-2, .2945E-3, .3492E-4, 

2 .3692E-5, .3524E-6,!. ,1. , .5, .1667, 

3 .4167E- 1,7*1. ,2. ,12. ,24. ,37.89, 

4 53.33,70.08,87.97,106.9,126.7, 

5 147.4,168.8,191.0,2.0,4.5,7.333, 

6 10 .42 , 13 .7, 7*1 ., 12 . 0, 24 .0, 37 . 89, 

7 53.33,70.08,87.97,106.9,126.7, 

8 147.4,168.8,191.0,1. ,3 .0,6.0, 

9 9.167,12.5,8*1./ 

FIRST EXECUTABLE STATEMENT 

GO TO (5,10), METH 
5 MAXDER = 12 

GOTO (15,20,25,30,35,40,45,50,55,60,65,70), NQ 
10 MAXDER = 5 



GOTO (75,80,85,90,95), NQ 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



THE FOLLOWING COEFFICIENTS SHOULD BE 
DEFINED TO MACHINE ACCURACY. FOR A 
GIVEN ORDER NQ, THEY CAN BE 
CALCULATED BY USE OF THE 
GENERATING POLYNOMIAL L (T) , WHOSE 
COEFFICIENTS ARE EL (I) .. L (T) = 
EL(1) + EL(2) *T + ... + 

EL(NQ+1) *T**NQ. FOR THE IMPLICIT 
ADAMS METHODS, L(T) IS GIVEN BY 
DL/DT = (T+1) * (T+2) * . . . 
*(T+NQ-1)/K, L(-l) =0, WHERE K = 



DGRCOOlO 

DGRC0020 

•DGRC0030 

DGRC0040 

DGRC0050 

DGRC0060 

DGRC0070 

DGRC0080 

DGRC0090 

DGRCOlOO 

DGRCOllO 

DGRC0120 

DGRC0130 

DGRC0140 

DGRC0150 

DGRC0160 

DGRC0170 

DGRC0180 

DGRC0190 

DGRC0200 

DGRC0210 

DGRC0220 

DGRC0230 

DGRC0240 

DGRC0250 

-DGRC0260 

DGRC0270 

DGRC0280 

DGRC0290 

DGRC0300 

DGRC0310 

DGRC0320 

DGRC0330 

DGRC0340 

DGRC0350 

DGRC0360 

DGRC0370 

DGRC0380 

DGRC0390 

DGRC0400 

DGRC0410 

DGRC0420 

DGRC0430 

DGRC0440 

DGRC0450 

DGRC0460 

DGRC0470 

DGRC0480 

DGRC0490 

DGRC0500 

DGRC0510 

DGRC0520 

DGRC0530 

DGRC0540 

DGRC0550 

DGRC0560 

DGRC0570 

DGRC0580 

DGRC0590 

DGRC0600 

DGRC0610 

DGRC0620 
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c 

c 

c 

c 

c 

c 

c 

c 

c 



FACTORIAL (NQ-1) . FOR THE GEAR DGRC063 0 
METHODS, L(T) = (T+1) * (T+2) * ... DGRC0640 
* (T+NQ) /K, WHERE K = DGRC0650 
FACTORIAL (NQ) * (1 + l/2 + ... + DGRC0660 
1/NQ) . THE ORDER IN WHICH THE DGRC0670 
GROUPS APPEAR BELOW IS.. IMPLICIT DGRC0680 
ADAMS METHODS OF ORDERS 1 TO 12 , DGRC0690 
BACKWARD DIFFERENTIATION METHODS DGRC0700 
OF ORDERS 1 TO 5. DGRC0710 



15 


EL(1) 


= l.ODO 


DGRC0720 




GO TO 


100 


DGRC0730 


20 


EL(1) 


= 0.5D0 


DGRC0740 




EL(3) 


= 0.5D0 


DGRC0750 




GO TO 


100 


DGRC0760 


25 


EL(1) 


= 4 .166666666666667D-01 


DGRC0770 




EL(3) 


= 0.75D0 


DGRC0780 




EL (4) 


= 1.666666666666667D-01 


DGRC0790 




GO TO 


100 


DGRC0800 


30 


EL(1) 


= 0.375D0 


DGRC0810 




EL(3) 


= 9 .166666666666667D-01 


DGRC0820 




EL (4) 


= 3 .333333333333333D-01 


DGRC0830 




EL(5) 


= 4 .166666666666667D-02 


DGRC0840 




GO TO 


100 


DGRC0850 


35 


EL(1) 


= 3 .486111111111111D-01 


DGRC0860 




EL(3) 


= 1.041666666666667D0 


DGRC0870 




EL (4) 


= 4 .861111111111111D-01 


DGRC0880 




EL(5) 


= 1.041666666666667D-01 


DGRC0890 




EL (6) 


= 8.333333333333333D-03 


DGRC0900 




GO TO 


100 


DGRC0910 


40 


EL(1) 


= 3 .298611111111111D-01 


DGRC0920 




EL(3) 


= 1.141666666666667D+00 


DGRC0930 




EL(4) 


= 0.625D+00 


DGRC0940 




EL(5) 


= 1. 770833333333333D-01 


DGRC0950 




EL(6) 


= 0.025D+00 


DGRC0960 




EL(7) 


= 1.388888888888889D-03 


DGRC0970 




GO TO 


100 


DGRC0980 


45 


EL(1) 


= 3 .155919312169312D-01 


DGRC0990 




EL(3) 


= 1.225D+00 


DGRCIOOO 




EL (4) 


= 7 .518518518518519D-01 


DGRClOlO 




EL(5) 


= 2 .552083333333333D-01 


DGRC1020 




EL(6) 


= 4 .861111111111111D-02 


DGRC1030 




EL(7) 


= 4 .861111111111111D-03 


DGRC1040 




EL (8) 


= 1 .984126984126984D-04 


DGRC1050 




GO TO 


100 


DGRC1060 


50 


EL(1) 


= 3 .042245370370370D-01 


DGRC1070 




EL(3) 


= 1 .296428571428571D+00 


DGRC1080 




EL (4) 


= 8 . 685185185185185D-01 


DGRC1090 




EL(5) 


= 3 .357638888888889D-01 


DGRCllOO 




EL(6) 


= 7.777777777777778D-02 


DGRClllO 




EL(7) 


= 1 .064814814814815D-02 


DGRC1120 




EL(8) 


= 7 .936507936507937D-04 


DGRC1130 




EL(9) 


= 2.480158730158730D-05 


DGRC1140 




GO TO 


100 


DGRC1150 


55 


EL(1) 


= 2 .948680004409171D-01 


DGRC1160 




EL(3) 


= 1 .358928571428571D+00 


DGRC1170 




EL(4) 


= 9 .765542328042328D-01 


DGRC1180 




EL(5) 


= 4 .171875D-01 


DGRC1190 




EL(6) 


= 1.113541666666667D-01 


DGRC1200 




EL(7) 


= 0.01875D+00 


DGRC1210 




EL(8) 


= 1 .934523809523810D-03 


DGRC1220 




EL(9) 


= 1.116071428571429D-04 


DGRC1230 




EL(10)= 2 .755731922398589D-06 


DGRC1240 



68 





GO TO 


100 


DGRC1250 


60 


EL(1) 


= 2 .869754464285714D-01 


DGRC1260 




EL(3) 


= 1 .414484126984127D+00 


DGRC1270 




EL(4) 


= 1 . 077215608465609D + 00 


DGRC1280 




EL(5) 


= 4 .985670194003527D-01 


DGRC1290 




EL (6) 


= 1.484375D-01 


DGRC1300 




EL(7) 


= 2 .906057098765432D-02 


DGRC1310 




EL (8) 


= 3 .720238095238095D-03 


DGRC1320 




EL(9) 


= 2 .996858465608466D-04 


DGRC1330 




EL(IO) 


= 1 .377865961199295D-05 


DGRC1340 




EL (11) 


= 2.755731922398589D-07 


DGRC1350 




GO TO 


100 


DGRC1360 


65 


EL(1) 


= 2 .801895964439367D-01 


DGRC1370 




EL(3) 


= 1.464484126984127D+00 


DGRC1380 




EL(4) 


= 1.171514550264550D+00 


DGRC1390 




EL(5) 


= 5.793581900352734D-01 


DGRC1400 




EL(6) 


= 1.883228615520282D-01 


DGRC1410 




EL(7) 


= 4 . 143036265432099D-02 


DGRC1420 




EL(8) 


= 6 .211144179894180D-03 


DGRC1430 




EL(9) 


= 6 .252066798941799D-04 


DGRC1440 




EL(IO) 


= 4.041740152851264D-05 


DGRC1450 




EL (11) 


= 1 .515652557319224D-06 


DGRC1460 




EL(12) 


= 2 .505210838544172D-08 


DGRC1470 




GO TO 


100 


DGRC1480 


70 


EL(1) 


= 2.742655400315991D-01 


DGRC1490 




EL(3) 


= 1 .509938672438672D+00 


DGRC1500 




EL (4) 


= 1 .260271164021164D+00 


DGRC1510 




EL(5) 


= 6 .592341820987654D-01 


DGRC1520 




EL(6) 


= 2 .304580026455027D-01 


DGRC1530 




EL(7) 


= 5 .569724610523222D-02 


DGRC1540 




EL (8) 


= 9 .439484126984127D-03 


DGRC1550 




EL (9) 


= 1.119274966931217D-03 


DGRC1560 




EL(IO) 


= 9 . 093915343915344D-05 


DGRC1570 




EL(ll) 


= 4 . 822530864197531D-06 


DGRC1580 




EL(12) 


= 1 .503126503126503D-07 


DGRC1590 




EL (13) 


= 2.087675698786810D-09 


DGRC1600 




GO TO 


100 


DGRC1610 

DGRC1620 


75 


EL(1) 


= l.OD+00 


DGRC1630 




GO TO 


100 


DGRC1640 


80 


EL(1) 


= 6 .666666666666667D-01 


DGRC1650 




EL(3) 


= 3 .333333333333333D-01 


DGRC1660 




GO TO 


100 


DGRC1670 


85 


EL(1) 


= 5 .454545454545455D-01 


DGRC1680 




EL(3) 


= EL(1) 


DGRC1690 




EL (4) 


= 9 .090909090909091D-02 


DGRC1700 




GO TO 


100 


DGRC1710 


90 


EL(1) 


= 0.48D+00 


DGRC1720 




EL(3) 


= 0.7D+00 


DGRC1730 




EL (4) 


= 0.2D+00 


DGRC1740 




EL(5) 


= 0.02D+00 


DGRC1750 




GO TO 


100 


DGRC1760 


95 


EL(1) 


= 4 .379562043795620D-01 


DGRC1770 




EL(3) 


= 8 . 211678832116788D-01 


DGRC1780 




EL(4) 


= 3 .102189781021898D-01 


DGRC1790 




EL(5) 


= 5 .474452554744526D-02 


DGRC1800 




EL (6) 


= 3 .649635036496350D-03 


DGRC1810 

DGRC1820 


100 


DO 105 K=l,3 


DGRC1830 




TQ(K) = PERTST (NQ,METH,K) 


DGRC1840 


105 


CONTINUE 


DGRC1850 




TQ(4) 


= .5D0*TQ(2) / (NQ+2) 


DGRC1860 
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RETURN 

END 



DGRC1870 

DGRC1880 
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IMSL ROUTINE NAME 



COMPUTER 
LATEST REVISION 
PURPOSE 

PRECISION/HARDWARE 



REQD. IMSL ROUTINES 
NOTATION 



COPYRIGHT 

WARRANTY 



SUBROUTINE DGRPS 



INTEGER 

DOUBLE PRECISION 



INTEGER 

★ 

★ 

REAL 

DOUBLE PRECISION 

★ 

COMMON /DBAND/ 
COMMON /GEAR/ 

★ 

★ 



- DGRPS 



- IBM/DOUBLE 

- NOVEMBER 1, 1984 

- NUCLEUS CALLED ONLY BY IMSL SUBROUTINE DGEAR 



- SINGLE AND DOUBLE /H3 2 

- SINGLE/H36,H48,H60 

- LUDATF,LEQT1B,UERTST,UGETI0 

- INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 

- 1984 BY IMSL, INC. ALL RIGHTS RESERVED. 

- IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 

APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 



( FCN , FCN J , Y , NO , CON , MITER , YMAX , SAVE 1 , SAVE 2 , PW , 
EQUIL, IPIV, lER) 

SPECIFICATIONS FOR ARGUMENTS 
NO, MITER, IPIV(l) , lER 

Y(N0, 1) ,CON, YMAX(l) , SAVEl (1) , SAVE2 (1) , PW(1) , 
EQUIL (1) 

SPECIFICATIONS FOR LOCAL VARIABLES 

NC , MFC , KFLAG , JSTART , NQUSED , NSTEP , NFE , NJE , NPW , 
NSQ , I , J1 , J , NERROR , NSAVEl , NSAVE2 , NEQUIL , NY , 
IDUMMY(23) ,NLIM, II , IJ, LIMl , LIM2 , NB , NLC , NUC , NWK 
SDUMMY(4) 

T , H , HMIN , HMAX , EPSC , UROUND , EPS J , HUSED , D , RO , YJ , R 
D1,D2,WA, DUMMY (40) 

NLC, NUC 

T , H , HMIN , HMAX , EPSC , UROUND , EPS J , HUSED , DUMMY , 
SDUMMY , NC , MFC , KFLAG , JSTART , NSQ , NQUSED , NSTEP , 
NFE , NJE , NPW, NERROR, NSAVEl , NSAVE2 , NEQUIL, NY, 
IDUMMY 

THIS ROUTINE IS CALLED BY DGRST TO 
COMPUTE AND PROCESS THE MATRIX P = 
I - H*EL(1)*J , WHERE J IS AN 
APPROXIMATION TO THE JACOBIAN. J 
IS COMPUTED, EITHER BY THE USER- 
SUPPLIED ROUTINE FCNJ IF MITER = 

1, OR BY FINITE DIFFERENCING IF 
MITER = 2 . J IS STORED IN PW AND 
REPLACED BY P, USING CON = 

-H*EL(1) . THEN P IS SUBJECTED TO 
LU DECOMPOSITION IN PREPARATION 
FOR LATER SOLUTION OF LINEAR 
SYSTEMS WITH P AS COEFFICIENT 
MATRIX. IN ADDITION TO VARIABLES 
DESCRIBED PREVIOUSLY, 
COMMUNICATION WITH DGRPS USES THE 



DGRPOOlO 

DGRP0020 

■DGRP0030 

DGRP0040 

DGRP0050 

DGRP0060 

DGRP0070 

DGRP0080 

DGRP0090 

DGRPOlOO 

DGRPOllO 

DGRP0120 

DGRP0130 

DGRP0140 

DGRP0150 

DGRP0160 

DGRP0170 

DGRP0180 

DGRP0190 

DGRP0200 

DGRP0210 

DGRP0220 

DGRP0230 

DGRP0240 

DGRP0250 

-DGRP0260 

DGRP0270 

DGRP0280 

DGRP0290 

DGRP0300 

DGRP0310 

DGRP0320 

DGRP0330 

DGRP0340 

DGRP0350 

DGRP0360 

DGRP0370 

DGRP0380 

DGRP0390 

,DGRP0400 

DGRP0410 

DGRP0420 

DGRP0430 

DGRP0440 

DGRP0450 

DGRP0460 

DGRP0470 

DGRP0480 

DGRP0490 

DGRP0500 

DGRP0510 

DGRP0520 

DGRP0530 

DGRP0540 

DGRP0550 

DGRP0560 

DGRP0570 

DGRP0580 

DGRP0590 

DGRP0600 

DGRP0610 

DGRP0620 
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FOLLOWING EPSJ = DSQRT (GROUND) , 
USED IN THE NUMERICAL JACOBIAN 
INCREMENTS . 

FIRST EXECUTABLE STATEMENT 

IF (NLC.EQ.-l) GO TO 45 

BANDED JACOBIAN CASE 

NB = NLC+NUC+1 

NWK = NB*N0+1 

IF (MITER. EQ. 2) GO TO 15 

MITER = 1 

NLIM = NB*N0 
DO 5 1=1, NLIM 
PW(I) = O.ODO 
5 CONTINUE 



CALL FCNJ(NC,T,Y,PW) 





DO 


10 


1=1 


,NLIM 






PW 


(I) 


= PW(I)*CON 


10 


CONTINUE 






GO 


TO 


35 




15 


D : 


= 0. 


. ODO 






DO 


20 


1=1 


,NC 


20 


D : 


= D+SAVE2 (I) **2 



RO = DABS (H) *DSQRT(D) *1 ,0D+03*UROUND 
DO 30 J=1,NC 
YJ = Y(J,1) 

R = EPSJ*YMAX(J) 

R = DMAXl (R,R0) 

Y(J,1) = Y(J,1)+R 
D = CON/R 

CALL FCN(NC,T,Y, SAVED 
LIMl = MAXO (1, J-NUC) 

LIM2 = MINO (NO, J+NLC) 

DO 25 I=LIM1,LIM2 

IJ = (J-I+NLC) *N0+I 
PW(IJ) = (SAVEl (I) -SAVE2 (I) ) *D 
25 CONTINUE 

Y(J,1) = YJ 
30 CONTINUE 

ADD IDENTITY MATRIX. 

35 DO 40 1=1, NC 

II = NLC*N0+I 
PW(II) = PW(II) +1.0D0 
40 CONTINUE 

DO LU DECOMPOSITION ON P 

CALL LEQT1B(PW,NC,NLC,NUC,N0,EQUIL,1,N0,1,PW(NWK) , lER) 
RETURN 

C FULL JACOBIAN CASE 

45 IF (MITER. EQ. 2) GO TO 55 
C MITER = 1 

CALL FCNJ(NC,T,Y,PW) 

DO 50 1=1, NSQ 
50 PW(I) = PW(I)*CON 
GO TO 75 

C MITER = 2 

55 D = O.ODO 

DO 60 1=1, NC 
60 D = D+SAVE2 (I) **2 

RO = DABS(H)*DSQRT(D)*1.0D+03*UROUND 
J1 = 0 



DGRP0630 

DGRP0640 

DGRP0650 

DGRP0660 

DGRP0670 

DGRP0680 

DGRP0690 

DGRP0700 

DGRP0710 

DGRP0720 

DGRP0730 

DGRP0740 

DGRP0750 

DGRP0760 

DGRP0770 

DGRP0780 

DGRP0790 

DGRP0800 

DGRP0810 

DGRP0820 

DGRP0830 

DGRP0840 

DGRP0850 

DGRP0860 

DGRP0870 

DGRP0880 

DGRP0890 

DGRP0900 

DGRP0910 

DGRP0920 

DGRP0930 

DGRP0940 

DGRP0950 

DGRP0960 

DGRP0970 

DGRP0980 

DGRP0990 

DGRPIOOO 

DGRPlOlO 

DGRP1020 

DGRP1030 

DGRP1040 

DGRP1050 

DGRP1060 

DGRP1070 

DGRP1080 

DGRP1090 

DGRPllOO 

DGRPlllO 

DGRP1120 

DGRP1130 

DGRP1140 

DGRP1150 

DGRP1160 

DGRP1170 

DGRP1180 

DGRP1190 

DGRP1200 

DGRP1210 

DGRP1220 

DGRP1230 

DGRP1240 
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DO 70 J=1,NC 
YJ = Y(J,1) 

R = EPSJ*YMAX(J) 

R = DMAXl (R, RO) 

Y(J,1) = Y(J,1)+R 

D = CON/R 

CALL FCN(NC,T, Y,SAVE1) 

DO 65 1=1, NC 

65 PW(I+J1) = (SAVEl (I) -SAVE2 (I) ) *D 
Y(J,1) = YJ 

J1 = Jl+NO 
70 CONTINUE 

C ADD IDENTITY MATRIX. 

75 J = 1 

DO 80 1=1, NC 

PW(J) = PW(J)+1.0D0 
J = J+(N0+1) 

80 CONTINUE 

DO LU DECOMPOSITION ON P, 

CALL LUDATF (PW, PW, NC , NO , 0 , D1 , D2 , IPIV, EQUIL, WA, lER) 
RETURN 
END 



DGRP1250 

DGRP1260 

DGRP1270 

DGRP1280 

DGRP1290 

DGRP1300 

DGRP1310 

DGRP1320 

DGRP1330 

DGRP1340 

DGRP1350 

DGRP1360 

DGRP1370 

DGRP1380 

DGRP1390 

DGRP1400 

DGRP1410 

DGRP1420 

DGRP1430 

DGRP1440 

DGRP1450 

DGRP1460 

DGRP1470 
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c 

c 

c« 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 



IMSL ROUTINE NAME 



COMPUTER 
LATEST REVISION 
PURPOSE 

PRECISION/HARDWARE 

REQD. IMSL ROUTINES 
NOTATION 

COPYRIGHT 

WARRANTY 



SUBROUTINE DGRIN 
INTEGER 

DOUBLE PRECISION 
INTEGER 

1 

2 

REAL 

DOUBLE PRECISION 

1 

COMMON /GEAR/ 

1 

2 

3 



DO 5 I = 
YO (I) 
5 CONTINUE 



1,NC 
= Y(I,1) 



L = JSTART + 1 
S = (TOUT - T) /H 
SI = l.ODO 
DO 15 J = 2,L 
*S1 = S1*S 



- DGRIN 



- IBM/DOUBLE 

- JANUARY 1, 1978 

NUCLEUS CALLED ONLY BY IMSL SUBROUTINE DGEAR 

- SINGLE AND DOUBLE /H3 2 

- SINGLE/H36,H48,H60 

- NONE REQUIRED 

- INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 

- 19 78 BY IMSL, INC. ALL RIGHTS RESERVED. 



IMSL WARRANTS ONLY THAT 
APPLIED TO THIS CODE. 
EXPRESSED OR IMPLIED, 



DGRIOOlO 
DGRI0020 
•DGRI0030 
DGRI0040 
DGRI0050 
DGRI0060 
DGRI0070 
DGRI0080 
DGRI0090 
DGRIOlOO 
DGRIOllO 
DGRI0120 
DGRI0130 
DGRI0140 
DGRI0150 
DGRI0160 
DGRI0170 
DGRI0180 
DGRI0190 
DGRI0200 
DGRI0210 

IMSL TESTING HAS BEEN DGRI0220 



NO OTHER WARRANTY, 
IS APPLICABLE. 



(TOUT, Y,N0, YO) 

SPECIFICATIONS FOR ARGUMENTS 
NO 

TOUT, YO (NO) , Y(N0, 1) 

SPECIFICATIONS FOR LOCAL VARIABLES 
NC , MFC , KFLAG , I , L , J , JSTART , NSQ , NQUSED , NSTEP , 
NFE , NJE , NPW, NERROR, NSAVEl , NSAVE2 , NEQUIL, NY, 
IDUMMY(23) 

SDUMMY(4) 

T , H , HMIN , HMAX , EPSC , UROUND , EPS J , HUSED , S , SI , 
DUMMY (40) 

T , H , HMIN , HMAX , EPSC , UROUND , EPS J , HUSED , DUMMY , 
SDUMMY , NC , MFC , KFLAG , JSTART , NSQ , NQUSED , NSTEP , 
NFE , NJE , NPW, NERROR, NSAVEl , NSAVE2 , NEQUIL, NY, 

I DUMMY 

FIRST EXECUTABLE STATEMENT 



THIS SUBROUTINE COMPUTES 

VALUES OF THE DEPENDENT VARIABLE 
Y AND STORES THEM IN YO . THE 
INTERPOLATION IS TO THE 
POINT T = TOUT, AND USES THE 
NORDSIECK HISTORY ARRAY Y, AS 
FOLLOWS . . 

NQ 

Y0(I) = SUM Y(I, J+1) *S**J , 

J=0 

WHERE S = -(T-TOUD/H. 



DGRI0230 
DGRI0240 
DGRI0250 
■DGRI0260 
DGRI0270 
DGRI0280 
DGRI0290 
DGRI0300 
DGRI0310 
DGRI0320 
DGRI0330 
DGRI0340 
DGRI0350 
DGRI0360 
DGRI0370 
DGRI0380 
DGRI0390 
DGRI0400 
DGRI0410 
DGRI0420 
DGRI0430 
DGRI0440 
DGRI0450 
DGRI0460 
INTERPOLATEDDGRI 04 7 0 



DGRI0480 

DGRI0490 

DGRI0500 

DGRI0510 

DGRI0520 

DGRI0530 

DGRI0540 

DGRI0550 

DGRI0560 

DGRI0570 

DGRI0580 

DGRI0590 

DGRI0600 

DGRI0610 

DGRI0620 
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DO 10 I = 1,NC DGRI0630 

Y0(I) = Y0(I) + S1*Y(I,J) DGRI0640 

10 CONTINUE DGRI0650 

15 CONTINUE DGRI0660 

RETURN DGRI0G70 

END DGRI0680 



75 



c 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



IMSL ROUTINE NAME 



COMPUTER 
LATEST REVISION 
PURPOSE 

USAGE 

ARGUMENTS A 
LU 



N 

lA 

IDGT 



D1 

D2 

I PVT 
EQUIL 

WA 

lER 



- LUDATF 



- IBM/DOUBLE 

- JANUARY 1, 1978 

- L-U DECOMPOSITION BY THE CROUT ALGORITHM 

WITH OPTIONAL ACCURACY TEST. 

- CALL LUDATF (A, LU, N, lA, IDGT, D1 , D2 , IPVT, 

EQUIL, WA, lER) 

- INPUT MATRIX OF DIMENSION N BY N CONTAINING 

THE MATRIX TO BE DECOMPOSED. 

- REAL OUTPUT MATRIX OF DIMENSION N BY N 

CONTAINING THE L-U DECOMPOSITION OF A 
ROWWISE PERMUTATION OF THE INPUT MATRIX. 
FOR A DESCRIPTION OF THE FORMAT OF LU, SEE 
EXAMPLE . 

- INPUT SCALAR CONTAINING THE ORDER OF THE 

MATRIX A. 

- INPUT SCALAR CONTAINING THE ROW DIMENSION OF 

MATRICES A AND LU EXACTLY AS SPECIFIED IN 
THE CALLING PROGRAM. 

- INPUT OPTION. 

IF IDGT IS GREATER THAN ZERO, THE NON- ZERO 
ELEMENTS OF A ARE ASSUMED TO BE CORRECT TO 
IDGT DECIMAL PLACES. LUDATF PERFORMS AN 
ACCURACY TEST TO DETERMINE IF THE COMPUTED 
DECOMPOSITION IS THE EXACT DECOMPOSITION 
OF A MATRIX WHICH DIFFERS FROM THE GIVEN 
ONE BY LESS THAN ITS UNCERTAINTY. 

IF IDGT IS EQUAL TO ZERO, THE ACCURACY TEST 
IS BYPASSED. 

- OUTPUT SCALAR CONTAINING ONE OF THE TWO 

COMPONENTS OF THE DETERMINANT. SEE 
DESCRIPTION OF PARAMETER D2 , BELOW. 

- OUTPUT SCALAR CONTAINING ONE OF THE 

TWO COMPONENTS OF THE DETERMINANT. THE 
DETERMINANT MAY BE EVALUATED AS (Dl) (2*^D2) 

- OUTPUT VECTOR OF LENGTH N CONTAINING THE 

PERMUTATION INDICES. SEE DOCUMENT 
(ALGORITHM) . 

- OUTPUT VECTOR OF LENGTH N CONTAINING 

RECIPROCALS OF THE ABSOLUTE VALUES OF 
THE LARGEST (IN ABSOLUTE VALUE) ELEMENT 
IN EACH ROW. 

- ACCURACY TEST PARAMETER, OUTPUT ONLY IF 

IDGT IS GREATER THAN ZERO. 

SEE ELEMENT DOCUMENTATION FOR DETAILS. 

- ERROR PARAMETER. (OUTPUT) 

TERMINAL ERROR 

lER = 129 INDICATES THAT MATRIX A IS 
ALGORITHMICALLY SINGULAR. (SEE THE 
CHAPTER L PRELUDE) . 

WARNING ERROR 

lER =34 INDICATES THAT THE ACCURACY TEST 
FAILED. THE COMPUTED SOLUTION MAY BE IN 
ERROR BY MORE THAN CAN BE ACCOUNTED FOR 
BY THE UNCERTAINTY OF THE DATA. THIS 



LUDAOOlO 
LUDA0020 
-LUDA0030 
LUDA0040 
LUDA0050 
LUDA0060 
LUDA0070 
LUDA0080 
LUDA0090 
LUDAOlOO 
LUDAOllO 
LUDA0120 
LUDA0130 
LUDA0140 
LUDA0150 
LUDA0160 
LUDA0170 
LUDA0180 
LUDA0190 
LUDA0200 
LUDA0210 
LUDA0220 
LUDA0230 
LUDA024 0 
LUDA0250 
LUDA0260 
LUDA0270 
LUDA0280 
LUDA0290 
LUDA0300 
LUDA0310 
LUDA0320 
LUDA0330 
LUDA0340 
LUDA0350 
LUDA0360 
LUDA0370 
LUDA0380 
LUDA0390 
LUDA0400 
LUDA0410 
.LUDA0420 
LUDA0430 
LUDA0440 
LUDA0450 
LUDA0460 
LUDA0470 
LUDA0480 
LUDA0490 
LUDA0500 
LUDA0510 
LUDA0520 
LUDA0530 
LUDA0540 
LUDA0550 
LUDA0560 
LUDA0570 
LUDA0580 
LUDA0590 
LUDA0600 
LUDA0610 
LUDA0620 
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on o noooooooooooonoonnnnooonn 



PRECISION/HARDWARE 



WARNING CAN BE PRODUCED ONLY IF IDGT IS 
GREATER THAN 0 ON INPUT. SEE CHAPTER L 
PRELUDE FOR FURTHER DISCUSSION. 

SINGLE AND DOUBLE /H3 2 
S INGLE /H3 6 , H4 8 , H60 



REQD. IMSL ROUTINES - UERTST,UGETIO 

NOTATION - INFORJMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 

REMARKS A TEST FOR SINGULARITY IS MADE AT TWO LEVELS: 

1. A ROW OF THE ORIGINAL MATRIX A IS NULL. 

2. A COLUMN BECOMES NULL IN THE FACTORIZATION PROCESS 

COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. 

WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 

APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 



SUBROUTINE LUDATF (A, LU, N, lA, IDGT, D1 , D2 , IPVT, EQUIL, WA, lER) 

DIMENSION A(IA,1) ,LU(IA,1) ,IPVT(1) , EQUIL (1) 

DOUBLE PRECISION A, LU, D1 , D2 , EQUIL, WA, ZERO, ONE , FOUR, SIXTN, SIXTH, 

* RN,WREL,BIGA, BIG, P, SUM, AI,WI,T, TEST, Q 

DATA ZERO, ONE, FOUR, SIXTN, SIXTH/0. DO, 1. DO, 4 .DO, 

* 16 .DO, .0625D0/ 

FIRST EXECUTABLE STATEMENT 
INITIALIZATION 

lER = 0 
RN = N 
WREL = ZERO 
D1 = ONE 
D2 = ZERO 
BIGA = ZERO 
DO 10 1=1, N 
BIG = ZERO 
DO 5 J=1,N 
P = A(I,J) 

LU(I,J) = P 
P = DABS(P) 

IF (P .GT. BIG) BIG = P 
5 CONTINUE 

IF (BIG .GT. BIGA) BIGA = BIG 
IF (BIG .EQ. ZERO) GO TO 110 
EQUIL (I) = ONE /BIG 
10 CONTINUE 

DO 105 J=1,N 
JMl = J-1 

IF (JMl .LT. 1) GO TO 40 

COMPUTE U( I, J) , 1 = 1 ,..., J-1 

DO 35 1 = 1, JMl 
SUM = LU(I,J) 

IMl = I-l 

IF (IDGT .EQ. 0) GO TO 25 

WITH ACCURACY TEST 

AI = DABS (SUM) 



LUDA0630 

LUDA0640 

LUDA0650 

LUDA0660 

LUDA0670 

LUDA0680 

LUDA0690 

LUDA0700 

LUDA0710 

LUDA0720 

LUDA0730 

LUDA0740 

LUDA0750 

LUDA0760 

LUDA0770 

.LUDA0780 

LUDA0790 

LUDA0800 

LUDA0810 

LUDA0820 

LUDA0830 

LUDA0840 

LUDA0850 

-LUDA0860 

LUDA0870 

LUDA0880 

LUDA0890 

LUDA0900 

LUDA0910 

LUDA0920 

LUDA0930 

LUDA0940 

LUDA0950 

LUDA0960 

LUDA0970 

LUDA0980 

LUDA0990 

LUDAIOOO 

LUDAlOlO 

LUDA1020 

LUDA1030 

LUDA1040 

LUDA1050 

LUDA1060 

LUDA1070 

LUDA1080 

LUDA1090 

LUDAllOO 

LUDAlllO 

LUDA1120 

LUDA1130 

LUDA1140 

LUDA1150 

LUDA1160 

LUDA1170 

LUDA1180 

LUDA1190 

LUDA1200 

LUDA1210 

LUDA1220 

LUDA1230 

LUDA1240 
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WI = ZERO 




LUDA1250 




IF (IMl .LT. 1) GO TO 20 




LUDA1260 




DO 15 K=1,IM1 




LUDA1270 




T = LU(I,K) *LU(K, J) 




LUDA1280 




SUM = SUM-T 




LUDA1290 




WI = WI+DABS(T) 




LUDA1300 


15 


CONTINUE 




LUDA1310 




LU(I,J) = SUM 




LUDA1320 


20 


WI = WI+DABS(SUM) 




LUDA1330 




IF (AI .EQ. ZERO) AI = BIGA 




LUDA1340 




TEST = WI/AI 




LUDA1350 




IF (TEST .GT. WREL) WREL = TEST 




LUDA1360 




GO TO 35 




LUDA1370 




WITHOUT ACCURACY 


LUDA1380 


25 


IF (IMl .LT. 1) GO TO 35 




LUDA1390 




DO 30 K=1,IM1 




LUDA1400 




SUM = SUM-LU(I,K) *LU(K, J) 




LUDA1410 


30 


CONTINUE 




LUDA1420 




LU(I,J) = SUM 




LUDA1430 


35 


CONTINUE 




LUDA1440 


40 


P = ZERO 




LUDA1450 




COMPUTE U(J,J) 


AND L(I, J) , I=J+1, . . 


. ,LUDA1460 




DO 70 I=J,N 




LUDA1470 




SUM = LU(I,J) 




LUDA1480 




IF (IDGT .EQ. 0) GO TO 55 




LUDA1490 




WITH ACCURACY 


TEST 


LUDA1500 




AI = DABS (SUM) 




LUDA1510 




WI = ZERO 




LUDA1520 




IF (JMl .LT. 1) GO TO 50 




LUDA1530 




DO 45 K=1,JM1 




LUDA1540 




T = LU(I,K) *LU(K, J) 




LUDA1550 




SUM = SUM-T 




LUDA1560 




WI = WI+DABS(T) 




LUDA1570 


45 


CONTINUE 




LUDA1580 




LU(I,J) = SUM 




LUDA1590 


50 


WI = WI+DABS(SUM) 




LUDA1600 




IF (AI .EQ. ZERO) AI = BIGA 




LUDA1610 




TEST = WI/AI 




LUDA1620 




IF (TEST .GT. WREL) WREL = TEST 




LUDA1630 




GO TO 65 




LUDA1640 




WITHOUT ACCURACY TEST 


LUDA1650 


55 


IF (JMl .LT. 1) GO TO 65 




LUDA1660 




DO 60 K=1,JM1 




LUDA1670 




SUM = SUM-LU(I,K) *LU(K, J) 




LUDA1680 


60 


CONTINUE 




LUDA1690 




LU(I,J) = SUM 




LUDA1700 


65 


Q = EQUIL (I) *DABS (SUM) 




LUDA1710 




IF (P .GE. Q) GO TO 70 




LUDA1720 




P = Q 




LUDA1730 




IMAX = I 




LUDA1740 


70 


CONTINUE 




LUDA1750 




TEST FOR ALGORITHMIC SINGULARITY 


LUDA1760 




IF (RN+P .EQ. RN) GO TO 110 




LUDA1770 




IF (J .EQ. IMAX) GO TO 80 




LUDA1780 




INTERCHANGE ROWS J AND IMAX 


LUDA1790 




D1 = -D1 




LUDA1800 




DO 75 K=1,N 




LUDA1810 




P = LU(IMAX,K) 




LUDA1820 




LU(IMAX,K) = LU(J,K) 




LUDA1830 




LU(J,K) = P 




LUDA1840 


75 


CONTINUE 




LUDA1850 




‘EQUIL (IMAX) = EQUIL (J) 




LUDA1860 



78 



80 IPVT(J) = IMAX 
D1 = D1*LU(J,J) 

85 IF (DABS(Dl) .LE. ONE) GO TO 90 
D1 = D1*SIXTH 
D2 = D2+FOUR 
GO TO 85 

90 IF (DABS(Dl) .GE. SIXTH) GO TO 95 
D1 = D1*SIXTN 
D2 = D2-FOUR 
GO TO 90 
95 CONTINUE 
JPl = J+1 

IF (JPl .GT. N) GO TO 105 

C DIVIDE BY PIVOT ELEMENT U(J,J) 

P = LU(J,J) 

DO 100 I=JP1,N 

LU(I, J) = LU(I, J) /P 
100 CONTINUE 
105 CONTINUE 

C PERFORM ACCURACY TEST 

IF (IDGT .EQ. 0) GO TO 9005 
P = 3*N+3 
WA = P*WREL 

IF (WA+IO.DO** (-IDGT) .NE. WA) GOTO 9005 
lER = 34 
GO TO 9000 

C ALGORITHMIC SINGULARITY 

110 lER = 129 
D1 = ZERO 
D2 = ZERO 
9000 CONTINUE 

C PRINT ERROR 

CALL UERTST (lER, 6HLUDATF) 

9005 RETURN 
END 



LUDA1870 

LUDA1880 

LUDA1890 

•LUDA1900 

LUDA1910 

LUDA1920 

LUDA1930 

LUDA1940 

LUDA1950 

LUDA1960 

LUDA1970 

LUDA1980 

LUDA1990 

LUDA2000 

LUDA2010 

LUDA2020 

LUDA2030 

LUDA2040 

LUDA2050 

LUDA2060 

LUDA2070 

LUDA2080 

LUDA2090 

LUDA2100 

LUDA2110 

LUDA2120 

LUDA2130 

LUDA2140 

LUDA2150 

LUDA2160 

LUDA2170 

LUDA2180 

LUDA2190 

LUDA2200 

LUDA2210 
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o n 
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IMSL ROUTINE NAME 



COMPUTER 
LATEST REVISION 
PURPOSE 

USAGE 

ARGUMENTS A 



B 

I PVT 

N 

lA 

X 

PRECISION/HARDWARE 

REQD. IMSL ROUTINES 
NOTATION 

COPYRIGHT 

WARRANTY 



SUBROUTINE LUELMF 

DIMENSION 
DOUBLE PRECISION 



DO 5 1=1, N 
5 X{I) = B{I) 

IW = 0 
DO 20 1=1, N 

IP = IPVT(I) 
SUM = X(IP) 
X(IP) = X(I) 

IF (IW .EQ. 0) 
IMl = I-l 



LUELMF 



IBM/DOUBLE 

• JANUARY 1, 1978 

ELIMINATION PART OF SOLUTION OF AX=B 
(FULL STORAGE MODE) 

CALL LUELMF (A, B , IPVT, N, lA, X) 

A = LU (THE RESULT COMPUTED IN THE IMSL 
ROUTINE LUDATF) WHERE L IS A LOWER 
TRIANGULAR MATRIX WITH ONES ON THE MAIN 
DIAGONAL. U IS UPPER TRIANGULAR. L AND U 
ARE STORED AS A SINGLE MATRIX A AND THE 
UNIT DIAGONAL OF L IS NOT STORED. (INPUT) 

B IS A VECTOR OF LENGTH N ON THE RIGHT HAND 
SIDE OF THE EQUATION AX=B . (INPUT) 

THE PERMUTATION MATRIX RETURNED FROM THE 
IMSL ROUTINE LUDATF, STORED AS AN N LENGTH 
VECTOR. (INPUT) 

ORDER OF A AND NUMBER OF ROWS IN B. (INPUT) 

ROW DIMENSION OF A EXACTLY AS SPECIFIED IN 
THE DIMENSION STATEMENT IN THE CALLING 
PROGRAM. (INPUT) 

THE RESULT X. (OUTPUT) 

- SINGLE AND DOUBLE/H32 

SINGLE/H36,H48,H60 

- NONE REQUIRED 

INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 

1978 BY IMSL, INC. ALL RIGHTS RESERVED. 

IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 
APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 



(A,B,IPVT,N,IA,X) 

A(IA,1) ,B(1) ,IPVT(1) ,X(1) 

A,B,X,SUM 

FIRST EXECUTABLE STATEMENT 
SOLVE LY = B FOR Y 



GO TO 15 



LUEFOOlO 

LUEF0020 

-LUEF0030 

LUEF0040 

LUEF0050 

LUEF0060 

LUEF0070 

LUEF0080 

LUEF0090 

LUEFOlOO 

LUEFOllO 

LUEF0120 

LUEF0130 

LUEF0140 

LUEF0150 

LUEF0160 

LUEF0170 

LUEF0180 

LUEF0190 

LUEF0200 

LUEF0210 

LUEF0220 

LUEF0230 

LUEF0240 

LUEF0250 

LUEF0260 

LUEF0270 

LUEF0280 

LUEF0290 

LUEF0300 

LUEF0310 

LUEF0320 

LUEF0330 

LUEF0340 

LUEF0350 

LUEF0360 

LUEF0370 

LUEF0380 

LUEF0390 

LUEF0400 

LUEF0410 

LUEF0420 

LUEF0430 

LUEF0440 

LUEF0450 

-LUEF0460 

LUEF0470 

LUEF0480 

LUEF0490 

LUEF0500 

LUEF0510 

LUEF0520 

LUEF0530 

LUEF0540 

LUEF0550 

LUEF0560 

LUEF0570 

LUEF0580 

LUEF0590 

LUEF0600 

LUEF0610 

LUEF0620 
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DO 10 J=IW,IM1 

SUM = SUM-A(I, J) *X(J) 

10 CONTINUE 

GO TO 20 

15 IF (SUM .NE. O.DO) IW = I 

20 X(I) = SUM 

C SOLVE UX 

DO 30 IB=1,N 
I = N+l-IB 
IPl = I+l 
SUM = X(I) 

IF (IPl .GT. N) GO TO 30 
DO 25 J=IP1,N 

SUM = SUM-A(I, J) *X(J) 

25 CONTINUE 
30 X(I) = SUM/A(I,I) 

RETURN 

END 



LUEF0630 

LUEF0640 

LUEF0650 

LUEF0660 

LUEF0670 

LUEF0680 

= Y FOR X LUEF0690 

LUEF0700 

LUEF0710 

LUEF0720 

LUEF0730 

LUEF0740 

LUEF0750 

LUEF0760 

LUEF0770 

LUEF0780 

LUEF0790 

LUEF0800 
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C IMSL ROUTINE NAME 
C 

C--- - 

c 

C COMPUTER 
C 

C LATEST REVISION 
C 

C PURPOSE 

C 

C 

C USAGE 

C 

C 



C ARGUMENTS A 
C 

C N 

C 

C NLC 

C 

C NUC 

C 

C lA 

C 

C 

C B 

C 

C 

C 

C 

C M 

C 

C IB 

C 

C 

C IJOB 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

C XL 

C 

c 

c 

c 

C lER 



- LEQTIB 



- IBM/DOUBLE 

- JANUARY 1, 1978 

- LINEAR EQUATION SOLUTION - BAND STORAGE 

MODE - SPACE ECONOMIZER SOLUTION 

- CALL LEQTIB (A, N, NLC , NUC , lA, B , M, IB , IJOB , XL, 

lER) 

- INPUT/OUTPUT MATRIX OF DIMENSION N BY 

(NUC+NLC+1) . SEE PARAMETER IJOB. 

- ORDER OF MATRIX A AND THE NUMBER OF ROWS IN 

B. (INPUT) 

- NUMBER OF LOWER CODIAGONALS IN MATRIX A. 

(INPUT) 

- NUMBER OF UPPER CODIAGONALS IN MATRIX A. 

(INPUT) 

- ROW DIMENSION OF MATRIX A EXACTLY AS 

SPECIFIED IN THE DIMENSION STATEMENT IN THE 
CALLING PROGRAM. (INPUT) 

- INPUT/OUTPUT MATRIX OF DIMENSION N BY M. 

ON INPUT, B CONTAINS THE M RIGHT-HAND SIDES 
OF THE EQUATION AX = B . ON OUTPUT, THE 
SOLUTION MATRIX X REPLACES B. IF IJOB = 1, 

B IS NOT USED. 

- NUMBER OF RIGHT HAND SIDES (COLUMNS IN B) . 

(INPUT) 

- ROW DIMENSION OF MATRIX B EXACTLY AS 

SPECIFIED IN THE DIMENSION STATEMENT IN THE 
CALLING PROGRAM. (INPUT) 

- INPUT OPTION PARAMETER. IJOB = I IMPLIES WHEN 

1 = 0, FACTOR THE MATRIX A AND SOLVE THE 
EQUATION AX = B. ON INPUT, A CONTAINS THE 
COEFFICIENT MATRIX OF THE EQUATION AX = B, 
WHERE A IS ASSUMED TO BE AN N BY N BAND 
MATRIX. A IS STORED IN BAND STORAGE MODE 
AND THEREFORE HAS DIMENSION N BY 
(NLC+NUC+1) . ON OUTPUT, A IS REPLACED 
BY THE U MATRIX OF THE L-U DECOMPOSITION 
OF A ROWWISE PERMUTATION OF MATRIX A. U 
IS STORED IN BAND STORAGE MODE. 

1 = 1, FACTOR THE MATRIX A. A CONTAINS THE 
SAME INPUT/OUTPUT INFORMATION AS IF 
IJOB = 0. 

1 = 2, SOLVE THE EQUATION AX = B . THIS 
OPTION IMPLIES THAT LEQTIB HAS ALREADY 
BEEN CALLED USING IJOB = 0 OR 1 SO THAT 
THE MATRIX A HAS ALREADY BEEN FACTORED. 

IN THIS CASE, OUTPUT MATRICES A AND XL 
MUST HAVE BEEN SAVED FOR REUSE IN THE 
CALL TO LEQTIB. 

- WORK AREA OF DIMENSION N* (NLC+1) . THE FIRST 

NLC*N LOCATIONS OF XL CONTAIN COMPONENTS OF 
THE L MATRIX OF THE L-U DECOMPOSITION OF A 
ROWWISE PERMUTATION OF A. THE LAST N 
LOCATIONS CONTAIN THE PIVOT INDICES. 

- ERROR PARAMETER. (OUTPUT) 



LElBOOlO 
LE1B0020 
-LE1B0030 
LE1B0040 
LE1B0050 
LE1B0060 
LE1B0070 
LE1B0080 
LE1B0090 
LElBOlOO 
LElBOllO 
LE1B0120 
LE1B0130 
LE1B0140 
LE1B0150 
LE1B0160 
LE1B0170 
LE1B0180 
LE1B0190 
LE1B0200 
LE1B0210 
LE1B0220 
LE1B0230 
LE1B0240 
LE1B0250 
LE1B0260 
LE1B0270 
LE1B0280 
LE1B0290 
LE1B0300 
LE1B0310 
LE1B0320 
LE1B0330 
LE1B0340 
LE1B0350 
LE1B0360 
LE1B0370 
LE1B0380 
, LE1B0390 
LE1B0400 
LE1B0410 
LE1B0420 
LE1B0430 
LE1B0440 
LE1B0450 
LE1B0460 
LE1B0470 
LE1B0480 
LE1B0490 
LE1B0500 
LE1B0510 
LE1B0520 
LE1B0530 
LE1B0540 
LE1B0550 
LE1B0560 
LE1B0570 
LE1B0580 
LE1B0590 
LE1B0600 
LE1B0610 
LE1B0620 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 



PRECISION/HARDWARE 



TERMINAL ERROR 

lER = 129 INDICATES THAT MATRIX A IS 
ALGORITHMICALLY SINGULAR. (SEE THE 
CHAPTER L PRELUDE) . 

- SINGLE AND DOUBLE /H3 2 
- SINGLE/H36,H48,H60 



REQD. IMSL ROUTINES - UERTST, UGETIO 



NOTATION 

COPYRIGHT 

W7VRRANTY 



- INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 

- 1978 BY IMSL, INC. ALL RIGHTS RESERVED. 

- IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 

APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 



SUBROUTINE LEQTIB (A, N, NLC , NUC , lA, B , M, IB , IJOB, XL, lER) 



DIMENSION 
DOUBLE PRECISION 
DATA 



A(IA,1) ,XL(N,1) ,B(IB,1) 

A,XL,B, P,Q,ZERO,ONE,RN 
ZERO/0. DO/, ONE/1 .ODO/ 

FIRST EXECUTABLE STATEMENT 



lER = 0 
JBEG = NLC+1 
NLCl = JBEG 

IF (IJOB .EQ. 2) GO TO 80 
RN = N 



RESTRUCTURE THE MATRIX 

FIND RECIPROCAL OF THE LARGEST 

ABSOLUTE VALUE IN ROW I 



1 = 1 

NC = JBEG+NUC 
NN = NC 
JEND = NC 

IF (N .EQ. 1 .OR. NLC . EQ . 0) GO TO 25 
5 K = 1 
P = ZERO 

DO 10 J = JBEG, JEND 
A(I,K) = A(I,J) 

Q = DABS (A(I,K) ) 

IF (Q .GT. P) P = Q 
K = K+1 
10 CONTINUE 
IF (P .EQ. 

XL (I, NLCl) 

IF (K .GT. 

DO 15 J = K,NC 
A(I,J) = ZERO 
15 CONTINUE 
20 I = I+l 

JBEG = JBEG-1 

IF (JEND -JBEG . EQ . N) JEND = JEND-1 
IF (I .LE. NLC) GO TO 5 
JBEG = I 
NN = JEND 
25 JElfJD = N-NUC 



ZERO) GO TO 135 

= ONE/P 

NC) GO TO 20 



LE1B0630 

LE1B0640 

LE1B0650 

LE1B0660 

LE1B0670 

LE1B0680 

LE1B0690 

LE1B0700 

LE1B0710 

LE1B0720 

LE1B0730 

LE1B0740 

LE1B0750 

LE1B0760 

LE1B0770 

LE1B0780 

LE1B0790 

LE1B0800 

LE1B0810 

LE1B0820 

-LE1B0830 

LE1B0840 

LE1B0850 

LE1B0860 

LE1B0870 

LE1B0880 

LE1B0890 

LE1B0900 

LE1B0910 

LE1B0920 

LE1B0930 

LE1B0940 

LE1B0950 

LE1B0960 

LE1B0970 

LE1B0980 

LE1B0990 

LEIBIOOO 

LElBlOlO 

LE1B1020 

LE1B1030 

LE1B1040 

LE1B1050 

LE1B1060 

LE1B1070 

LE1B1080 

LE1B1090 

LEIBIIOO 

LElBlllO 

LE1B1120 

LE1B1130 

LE1B1140 

LE1B1150 

LE1B1160 

LE1B1170 

LE1B1180 

LE1B1190 

LE1B1200 

LE1B1210 

LE1B1220 

LE1B1230 

LE1B1240 
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DO 40 I = JBEG,N 




LE1B1250 




P = ZERO 




LE1B1260 




DO 30 J = 1,NN 




LE1B1270 




Q = DABS (A(I, J) ) 




LE1B1280 




IF (Q ,GT. P) P = Q 




LE1B1290 


30 


CONTINUE 




LE1B1300 




IF (P .EQ. ZERO) GO TO 135 




LE1B1310 




XL(I,NLC1) = ONE/P 




LE1B1320 




IF (I .EQ. JEND) GO TO 37 




LE1B1330 




IF (I .LT. JEND) GO TO 40 




LE1B1340 




K = NN+1 




LE1B1350 




DO 35 J = K,NC 




LE1B1360 




A(I,J) = ZERO 




LE1B1370 


35 


CONTINUE 




LE1B1380 


37 


NN = NN-1 




LE1B1390 


40 


CONTINUE 




LE1B1400 




L = NLC 




LE1B1410 




L-U DECOMPOSITION 


LE1B1420 




DO 75 K = 1,N 




LE1B1430 




P = DABS (A(K, 1) ) *XL(K,NLC1) 




LE1B1440 




I = K 




LE1B1450 




IF (L .LT. N) L = L+1 




LE1B1460 




K1 = K+1 




LE1B1470 




IF (K1 .GT. L) GO TO 50 




LE1B1480 




DO 45 J = K1,L 




LE1B1490 




Q = DABS (A(J, 1) ) *XL(J,NLC1) 




LE1B1500 




IF (Q .LE. P) GO TO 45 




LE1B1510 




P = Q 




LE1B1520 




I = J 




LE1B1530 


45 


CONTINUE 




LE1B1540 


50 


XL(I,NLC1) = XL(K,NLC1) 




LE1B1550 




XL(K,NLC1) = I 




LE1B1560 




SINGULARITY 


FOUND 


LE1B1570 




Q = RN+P 




LE1B1580 




IF (Q .EQ. RN) GO TO 135 




LE1B1590 




INTERCHANGE 


ROWS I AND K 


LE1B1600 




IF (K .EQ. I) GO TO 60 




LE1B1610 




DO 55 J = 1,NC 




LE1B1620 




P = A(K,J) 




LE1B1630 




A(K,J) = A(I,J) 




LE1B1640 




A(I,J) = P 




LE1B1650 


55 


CONTINUE 




LE1B1660 


60 


IF (K1 .GT. L) GO TO 75 




LE1B1670 




DO 70 I = K1,L 




LE1B1680 




P = A(I,1) /A(K,1) 




LE1B1690 




IK = I-K 




LE1B1700 




XL(K1,IK) = P 




LE1B1710 




DO 65 J = 2,NC 




LE1B1720 




A(I,J-1) = A(I, J) -P*A(K, J) 




LE1B1730 


65 


CONTINUE 




LE1B1740 




A (I, NO = ZERO 




LE1B1750 


70 


CONTINUE 




LE1B1760 


75 


CONTINUE 




LE1B1770 




IF (IJOB .EQ. 1) GO TO 9005 




LE1B1780 




FORWARD SUBSTITUTION 


LE1B1790 


80 


L = NLC 




LE1B1800 




DO 105 K = 1,N 




LE1B1810 




I = XL(K,NLC1) 




LE1B1820 




IF (I .EQ. K) GO TO 90 




LE1B1830 




DO 85 J = 1,M 




LE1B1840 




P = B(K,J) 




LE1B1850 




B(K,J) = B(I,J) 




LE1B1860 



84 



B(I,J) = P 
85 CONTINUE 

90 IF (L .LT. N) L = L+1 

K1 = K+1 

IF (K1 .GT. L) GO TO 105 
DO 100 I = K1,L 
IK = I-K 
P = XL(K1,IK) 

DO 95 J = 1,M 

B(I, J) = B(I, J) -P*B(K, J) 
95 CONTINUE 

100 CONTINUE 

105 CONTINUE 

C BACKW 

JBEG = NUC+NLC 
DO 125 J = 1,M 
L = 1 
K1 = N+1 
DO 120 I = 1,N 
K = Kl-I 
P = B(K,J) 

IF (L .EQ. 1) GO TO 115 
DO 110 KK = 2,L 
IK = KK+K 

P = P-A(K,KK) *B(IK-1, J) 
110 CONTINUE 

115 B(K,J) = P/A(K,1) 

IF (L .LE. JBEG) L = L+1 
120 CONTINUE 
125 CONTINUE 
GO TO 9005 
135 lER = 129 
9000 CONTINUE 

CALL UERTST (lER, 6HLEQT1B) 

9005 RETURN 
END 



LE1B1870 

LE1B1880 

LE1B1890 

LE1B1900 

LE1B1910 

LE1B1920 

LE1B1930 

LE1B1940 

LE1B1950 

LE1B1960 

LE1B1970 

LE1B1980 

LE1B1990 

SUBSTITUTION LE1B2000 

LE1B2010 

LE1B2020 

LE1B2030 

LE1B2040 

LE1B2050 

LE1B2060 

LE1B2070 

LE1B2080 

LE1B2090 

LE1B2100 

LE1B2110 

LE1B2120 

LE1B2130 

LE1B2140 

LE1B2150 

LE1B2160 

LE1B2170 

LE1B2180 

LE1B2190 

LE1B2200 

LE1B2210 

LE1B2220 
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IMSL ROUTINE NAME 



- UERTST 



COMPUTER 
LATEST REVISION 
PURPOSE 
USAGE 

ARGUMENTS lER 



NAME 

PRE C I S I ON /HARDWARE 
REQD. IMSL ROUTINES 
NOTATION 



- IBM/SINGLE 

- JUNE 1, 1982 

- PRINT A MESSAGE REFLECTING AN ERROR CONDITION 

- CALL UERTST (IER,NAME) 

- ERROR PARAMETER. (INPUT) 

lER = I+J WHERE 

I = 128 IMPLIES TERMINAL ERROR MESSAGE, 

I = 64 IMPLIES WARNING WITH FIX MESSAGE, 

I = 32 IMPLIES WARNING MESSAGE. 

J = ERROR CODE RELEVANT TO CALLING 
ROUTINE . 

- A CHARACTER STRING OF LENGTH SIX PROVIDING 

THE NAME OF THE CALLING ROUTINE. (INPUT) 



- SINGLE /ALL 

- UGETIO,USPKD 

INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 



REMARKS THE ERROR MESSAGE PRODUCED BY UERTST IS WRITTEN 

TO THE STANDARD OUTPUT UNIT. THE OUTPUT UNIT 
NUMBER CAN BE DETERMINED BY CALLING UGETIO AS 
FOLLOWS . . CALL UGETIO ( 1 , NIN, NOUT) . 

THE OUTPUT UNIT NUMBER CAN BE CHANGED BY CALLING 
UGETIO AS FOLLOWS . . 

NIN = 0 

NOUT = NEW OUTPUT UNIT NUMBER 
CALL UGETIO (3, NIN, NOUT) 

SEE THE UGETIO DOCUMENT FOR MORE DETAILS. 

COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. 

WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 

APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 



SUBROUTINE UERTST (IER,NAME) 

SPECIFICATIONS FOR ARGUMENTS 

lER 

NAME (1) 

SPECIFICATIONS FOR LOCAL VARIABLES 
I , lEQ , lEQDF , lOUNIT , LEVEL , LEVOLD , NAMEQ ( 6 ) , 
NAMSET(6) ,NAMUPK(6) ,NIN,NMTB 
NAMSET/IHU, IHE, IHR, IHS, IHE, IHT/ 

NAMEQ/6*1H / 

LEVEL^/, IEQDF/0/, IEQ/1H=/ 

UNPACK NAME INTO NAMUPK 
FIRST EXECUTABLE STATEMENT 
CALL USPKD (NAME, 6, NAMUPK, NMTB) 



INTEGER 

INTEGER 

INTEGER 

DATA 

DATA 

DATA 



UERTOOlO 

UERT0020 

■UERT0030 

UERT0040 

UERT0050 

UERT0060 

UERT0070 

UERT0080 

UERT0090 

UERTOlOO 

UERTOllO 

UERT0120 

UERT0130 

UERT0140 

UERT0150 

UERT0160 

UERT0170 

UERT0180 

UERT0190 

UERT0200 

UERT0210 

UERT0220 

UERT0230 

UERT0240 

UERT0250 

UERT0260 

UERT0270 

UERT0280 

UERT0290 

UERT0300 

UERT0310 

UERT0320 

UERT0330 

UERT0340 

UERT0350 

UERT0360 

UERT0370 

UERT0380 

UERT0390 

UERT0400 

UERT0410 

UERT0420 

UERT0430 

UERT0440 

UERT0450 

UERT0460 

UERT0470 

-UERT0480 

UERT0490 

UERT0500 

UERT0510 

UERT0520 

UERT0530 

UERT0540 

UERT0550 

UERT0560 

UERT0570 

UERT0580 

UERT0590 

UERT0600 

UERT0610 

UERT0620 
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o n o o 



GET OUTPUT UNIT NUMBER 



CHECK lER 



CALL UGETIO(l,NIN, lOUNIT) 

IF (IER.GT.999) GO TO 25 

IF (IER.LT.-32) GO TO 55 

IF (IER.LE.128) GO TO 5 

IF (LEVEL. LT.l) GO TO 30 

PRINT TERMINAL MESSAGE 

IF (lEQDF.EQ.l) WRITE ( lOUNIT, 35 ) lER, NAMEQ, lEQ, NAMUPK 
IF (lEQDF.EQ.O) WRITE ( lOUNIT, 35) IER,NAMUPK 
GO TO 30 

5 IF (IER.LE.64) GO TO 10 
IF (LEVEL.lt. 2) GO TO 30 

PRINT WARNING WITH FIX MESSAGE 
IF (lEQDF.EQ.l) WRITE ( lOUNIT , 4 0 ) lER, NAMEQ , lEQ , NAMUPK 
IF (lEQDF.EQ.O) WRITE (lOUNIT, 40 ) IER,NAMUPK 
GO TO 30 

10 IF (IER.LE.32) GO TO 15 

PRINT WARNING MESSAGE 

IF (LEVEL.lt. 3) GO TO 30 

IF (lEQDF.EQ.l) WRITE ( lOUNIT, 45 ) lER, NAMEQ, lEQ, NAMUPK 
IF (lEQDF.EQ.O) WRITE ( lOUNIT, 45 ) IER,NAMUPK 
GO TO 30 
15 CONTINUE 

CHECK FOR UERSET CALL 

DO 20 1=1,6 

IF (NAMUPK(I) .NE.NAMSET(I) ) GO TO 25 
20 CONTINUE 

LEVOLD = LEVEL 
LEVEL = lER 
lER = LEVOLD 

IF (LEVEL.lt. 0) LEVEL = 4 
IF (LEVEL.GT.4) LEVEL = 4 
GO TO 30 
25 CONTINUE 

IF (LEVEL.lt. 4) GO TO 30 

PRINT NON -DEFINED MESSAGE 
IF (lEQDF.EQ.l) WRITE ( lOUNIT, 50 ) lER, NAMEQ, lEQ, NAMUPK 
IF (lEQDF.EQ.O) WRITE ( lOUNIT, 50 ) lER, NAMUPK 
30 lEQDF = 0 
RETURN 

35 FORMAT (19H *** TERMINAL ERROR, 1 OX, 7H (lER = ,13, 

1 2 OH) FROM IMSL ROUTINE ,6A1,A1,6A1) 

40 FORMAT(27H *** WARNING WITH FIX ERROR, 2X , 7H ( lER = ,13, 

1 2 OH) FROM IMSL ROUTINE ,6A1,A1,6A1) 

45 FORMAT (18H WARNING ERROR, 1 IX, 7H (lER = , 13 , 

1 2 OH) FROM IMSL ROUTINE ,6A1,A1,6A1) 

50 FORMAT (20H UNDEFINED ERROR, 9X , 7H ( lER = ,15, 

1 20H) FROM IMSL ROUTINE ,6A1,A1,6A1) 



SAVE P FOR P = R CASE 
P IS THE PAGE NAMUPK 
R IS THE ROUTINE NAMUPK 

55 lEQDF = 1 
DO 60 1=1,6 

60 NAMEQ(I) = NAMUPK(I) 

65 RETURN 
END 



UERT0630 

UERT0640 

UERT0650 

■JERT0660 

UERT0670 

UERT0680 

UERT0690 

UERT0700 

UERT0710 

UERT0720 

UERT0730 

UERT0740 

UERT0750 

UERT0760 

UERT0770 

UERT0780 

UERT0790 

UERT0800 

UERT0810 

UERT0820 

UERT0830 

UERT0840 

UERT0850 

UERT0860 

UERT0870 

UERT0880 

UERT0890 

UERT0900 

UERT0910 

UERT0920 

UERT0930 

UERT0940 

UERT0950 

UERT0960 

UERT0970 

UERT0980 

UERT0990 

UERTIOOO 

UERTlOlO 

UERT1020 

UERT1030 

UERT1040 

UERT1050 

UERT1060 

UERT1070 

UERT1080 

UERT1090 

UERTllOO 

UERTlllO 

UERT1120 

UERT1130 

UERT1140 

UERT1150 

UERT1160 

UERT1170 

UERT1180 

UERT1190 

UERT1200 
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c 



IMSL ROUTINE NAME - UGETIO 



COMPUTER 
LATEST REVISION 
PURPOSE 



- IBM/ SINGLE 

- JUNE 1, 1981 

- TO RETRIEVE CURRENT VALUES AND TO SET NEW 
VALUES FOR INPUT AND OUTPUT UNIT 
IDENTIFIERS. 



USAGE 



CALL UGETIO (IOPT,NIN,NOUT) 



ARGUMENTS lOPT 



NIN 

NOUT 



OPTION PARAMETER. (INPUT) 

IF IOPT=l, THE CURRENT INPUT AND OUTPUT 
UNIT IDENTIFIER VALUES ARE RETURNED IN NIN 
AND NOUT, RESPECTIVELY. 

IF IOPT=2, THE INTERNAL VALUE OF NIN IS 
RESET FOR SUBSEQUENT USE. 

IF IOPT=3, THE INTERNAL VALUE OF NOUT IS 
RESET FOR SUBSEQUENT USE. 

INPUT UNIT IDENTIFIER. 

OUTPUT IF IOPT=l, INPUT IF IOPT=2 . 

OUTPUT UNIT IDENTIFIER. 

OUTPUT IF IOPT=l, INPUT IF IOPT=3 . 



PRECISION/HARDWARE - SINGLE/ALL 



REQD. IMSL ROUTINES - NONE REQUIRED 

NOTATION - INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 

REMARKS EACH IMSL ROUTINE THAT PERFORMS INPUT AND/OR OUTPUT 

OPERATIONS CALLS UGETIO TO OBTAIN THE CURRENT UNIT 
IDENTIFIER VALUES, IF UGETIO IS CALLED WITH IOPT=2 OR 
IOPT=3, NEW UNIT IDENTIFIER VALUES ARE ESTABLISHED. 
SUBSEQUENT INPUT/OUTPUT IS PERFORMED ON THE NEW UNITS. 



COPYRIGHT 



1978 BY IMSL, INC. ALL RIGHTS RESERVED. 



UGETOOlO 

UGET0020 

-UGET0030 

UGET0040 

UGET0050 

UGET0060 

UGET0070 

UGET0080 

UGET0090 

UGETOlOO 

UGETOllO 

UGET0120 

UGET0130 

UGET0140 

UGET0150 

UGET0160 

UGET0170 

UGET0180 

UGET0190 

UGET0200 

UGET0210 

UGET0220 

UGET0230 

UGET0240 

UGET0250 

UGET0260 

UGET0270 

UGET0280 

UGET0290 

UGET0300 

UGET0310 

UGET0320 

UGET0330 

UGET0340 

UGET0350 

UGET0360 

UGET0370 

UGET0380 

UGET0390 

UGET0400 

UGET0410 

UGET0420 

UGET0430 



WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN UGET044 0 

APPLIED TO THIS CODE. NO OTHER WARRANTY, UGET0450 

EXPRESSED OR IMPLIED, IS APPLICABLE. UGET0460 

UGET0470 



UGET0480 

UGET0490 

SUBROUTINE UGETIO ( lOPT, NIN, NOUT) UGET0500 

SPECIFICATIONS FOR ARGUMENTS UGET0510 



INTEGER 




lOPT, NIN, NOUT 


UGET0520 






SPECIFICATIONS FOR LOCAL VARIABLES 


UGET0530 


INTEGER 




NIND, NOUTD 


UGET0540 


DATA 




NIND/5/,NOUTD/6/ 


UGET0550 






FIRST EXECUTABLE STATEMENT 


UGET0560 


IF (IOPT.EQ.3) 


GO 


TO 10 


UGET0570 


IF (IOPT.EQ.2) 


GO 


TO 5 


UGET0580 


IF (lOPT.NE.l) 


GO 


TO 9005 


UGET0590 


NIN = NIND 






UGET0600 


NOUT = NOUTD 






UGET0610 


GO TO 9005 






UGET0620 
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5 NIND = NIN UGET0630 

GO TO 9005 UGET0640 

10 NOUTD = NOUT UGET0650 

9005 RETURN UGET0660 

Eiro UGET0670 
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IMSL ROUTINE NAME 



USPKD 



COMPUTER 
LATEST REVISION 
PURPOSE 



USAGE 

ARGUMENTS 



C 

C 

C 



- IBM/ SINGLE 

- NOVEMBER 1, 1984 

- NUCLEUS CALLED BY IMSL ROUTINES THAT HAVE 

CHARACTER STRING ARGUMENTS 

- CALL USPKD ( PACKED, NCHARS,UNPAKD,NCHMTB) 



PACKED - CHARACTER STRING TO BE UNPACKED. (INPUT) 
NCHARS - LENGTH OF PACKED. (INPUT) SEE REMARKS. 
UNPAKD - INTEGER ARRAY TO RECEIVE THE UNPACKED 
REPRESENTATION OF THE STRING. (OUTPUT) 
NCHMTB - NCHARS MINUS TRAILING BLANKS. (OUTPUT) 



PRECISION/HARDWARE - SINGLE /ALL 
REQD. IMSL ROUTINES - NONE 
REMARKS 1 . 

2 . 



USPKD UNPACKS A CHARACTER STRING INTO AN INTEGER ARRAY 
IN (Al) FORMAT. 

UP TO 129 CHARACTERS MAY BE USED. ANY IN EXCESS OF 
THAT ARE IGNORED. 



COPYRIGHT - 1984 BY IMSL, INC. ALh RIGHTS RESERVED. 

WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 

APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 



SUBROUTINE USPKD ( PACKED , NCHARS , UNPAKD , NCHMTB) 

SPECIFICATIONS FOR ARGUMENTS 
INTEGER NC , NCHARS , NCHMTB 

LOGICAL^ 1 UNPAKD ( 1 ) , PACKED ( 1 ) , LBYTE , LBLANK 

INTEGER*2 IBYTE , IBLANK 

EQUIVALENCE ( LBYTE , IBYTE ) 

DATA LBLANK /IH / 

DATA IBYTE /IH / 

DATA IBLANK /IH / 

INITIALIZE NCHMTB 



NCHMTB = 0 

IF (NCHARS. LE.O) RETURN 

NC = MINO (129, NCHARS) 
NWORDS = NC*4 
J = 1 

DO 110 I = 1, NWORDS, 4 
UNPAKD (I) = PACKED (J) 
UNPAKD (I+l) = LBLANK 
UNPAKD (1+2) = LBLANK 
UNPAKD (1+3) = LBLANK 
110 J = J+1 



DO* 200 N = 1, NWORDS, 4 



RETURN IF NCHARS IS LE ZERO 



SET NC=NUMBER OF CHARS TO BE DECODED 



CHECK UNPAKD ARRAY AND SET NCHMTB 
BASED ON TRAILING BLANKS FOUND 



USPKOOlO 

USPK0020 

■USPK0030 

USPK0040 

USPK0050 

USPK0060 

USPK0070 

USPK0080 

USPK0090 

USPKOlOO 

USPKOllO 

USPK0120 

USPK0130 

USPK0140 

USPK0150 

USPK0160 

USPK0170 

USPK0180 

USPK0190 

USPK0200 

USPK0210 

USPK0220 

USPK0230 

USPK0240 

USPK0250 

USPK0260 

USPK0270 

USPK0280 

USPK0290 

USPK0300 

USPK0310 

USPK0320 

USPK0330 

USPK0340 

-USPK0350 

USPK0360 

USPK0370 

USPK0380 

USPK0390 

USPK0400 

USPK0410 

USPK0420 

USPK0430 

USPK0440 

USPK0450 

USPK0460 

USPK0470 

USPK0480 

USPK0490 

USPK0500 

USPK0510 

USPK0520 

USPK0530 

USPK0540 

USPK0550 

USPK0560 

USPK0570 

USPK0580 

USPK0590 

USPK0600 

USPK0610 

USPK0620 
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NN = NWORDS - N - 2 
LBYTE = UNPAKD(NN) 

IFdBYTE .NE. IBLANK) GO TO 210 
200 CONTINUE 
NN = 0 

210 NCHMTB = (NN + 3) / 4 

RETURN 
END 



USPK0630 

USPK0640 

USPK0650 

USPK0660 

USPK0670 

USPK0680 

USPK0690 

USPK0700 
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